Online Ressource: Varroa Host Switch Genomics and Demography

This online ressource is there to provide additional interactive visualization tool for the Varroa mites sampling and sequencing data obtained from the associated paper “XXXX”. A Snakemake pipeline created and used for genomics mapping, variant calling and analysis is provided in this host GitHub repository. Most genomics and demographic analysis using fastsimcoal2 were computed on the Okinawa Institute of Science and Technology OIST cluster, Sango.

Click on “code” button on the right to make appear the libraries used and R code for each section.

library("tidyverse")
library("knitr") # for the markdown
library("kableExtra") # for creating a scrolling table
library("ggplot2") # for plotting 
library("maps")
library("leaflet")
library("htmltools")
library("rgdal")
library("gridExtra")
library("gghighlight")
library("adegenet")
library("vcfR")
library("ggmap")
library("poppr")
library("RColorBrewer")
library("gridExtra")
library("igraph")
library("StAMPP")
library("lattice")
library("reshape2")
library("ggrepel") # for text labels
library("pcadapt")
library("data.table")
library("ggthemes")
library("ggpubr")
library("boot")
library("plotly")
library("forcats")
library("plotly")

1. Distribution maps of the Varroa mites sampled

Here, we have 43 samples of Varroa mites for which we have sequenced the whole genome succesfully. We use the leaflet package to create an interactive map showing the location and from which host, each Varroa mite was collected.

ID = Name code of the samples obtained from the CSIRO varroa collection
Species = Identified species according to mtDNA of the mtDNA COX1 partial gene
Host collected = Varroa mites were collected from within honey bee colonies identified by the collector as either Western honey bee Apis mellifera or Eastern honey bee Apis cerana.
Location/Country/Date = Exact informations on sampling
Site = Code name of each sampling site corresponding on the map Approximate.X/Y = GPS coordinates from samples obtained before 2008 are infered from approximate location given. All samples geo-coordinates after 2008 were obtained directly from field collection.

COI_haplogroup = mtDNA haplogroup identified from mtDNA COX1 partial region (reconstruct from HiSeq4000 reads and confirmed by Sanger sequencing).
Haplotype = mtDNA haplotype name by concatenating four genes (COX1, COX3, ATP6 and ND2)

Host_read = Host DNA identified by number of reads mapped against either A. cerana or A. mellifera
reads_mellifera/cerana = Number of reads mapped against reference genome of A. cerana or A. mellifera (Q > 20)

# Import the metadata for the Varroa samples 
metadata <- read.csv("metadata_varroa_asia_15052020.csv", header = TRUE)
#head(metadata)

## Create a scrolling table
kable(cbind(metadata)) %>%
  kable_styling(bootstrap_options = "striped", font_size = 9)  %>%
  scroll_box(width = "100%", height = "400px") #%>%

ID	SEQ_ID	Species	Host_collected	Location	Country	Date	Site	Approximate.X	Approximate.Y	COI_haplogroup	Haplotype	Host_read	reads_cerana	reads_mellifera	Mean_Read_Depth	Total_reads	Mapped	Mapping_rate	Mapped_in_pair	Platform
V212	VD212	V. destructor	A. cerana	Gwangju, Jeonlabug Province	South Korea	01/1996	KOR02	35.159540	126.85260	VD Korea K1	VD Korea K1-1/K1-2	cerana	56	0	16.6	66272622	41762260	63.0	38388417	HiSeq 4000
V642_1	VD642_1	V. destructor	A. cerana	Zhuhai county, Guangdong Province	China	11/2001	CHN08	22.614010	112.63183	VD China C1	VD China C1-3	cerana	2857	7	47.1	119107269	116682672	98.0	112230978	HiSeq 4000
V657_1	VD657_1	V. destructor	A. cerana	Zhongshan, Guangdong Province	China	04/2002	CHN10	23.946090	114.69726	VD China C1	VD China C1-3	cerana	46	38	21.5	54511467	53363635	97.9	49849866	HiSeq 4000
V651_1	VD651_1	V. destructor	A. cerana	Dayao county, Yunnan Province	China	04/2002	CHN05	25.729510	101.33661	VD China C3	VD China C3-2	cerana	167	6	28.2	71261265	69804599	98.0	66308649	HiSeq 4000
V577_1	VD577_1	V. destructor	A. cerana	Yunnan Agricultural University, Kunming, Yunnan Province	China	05/2002	CHN06	25.125650	102.74318	VD China C4	VD China C4-2	cerana	125	4	35.0	91721740	89577965	97.7	85757545	HiSeq 4000
V654_1	VD654_1	V. destructor	A. cerana	Xishuangbanna, Yunnan Province	China	04/2002	CHN07	22.008810	100.79714	VD Viet Nam V1	VD Viet Nam V1-5	cerana	626	0	25.1	63453379	62161106	98.0	58776768	HiSeq 4000
V658_2	VD658_2	V. destructor	A. cerana	Xishuangbanna, Yunnan Province	China	04/2002	CHN07	22.008810	100.79714	VD Viet Nam V1	VD Viet Nam V1-6	cerana	2693	8	15.0	44339428	38654491	87.2	35187264	HiSeq 4000
V676_2	VD676_2	V. destructor	A. cerana	Chiang Mai University	Thailand	01/2003	THA01	18.805960	98.95336	VD Viet Nam V1	VD Viet Nam V1-7	cerana	366	10	5.9	29358273	19254785	65.6	17719549	HiSeq 4000
V149	VD149	V. destructor	A. cerana	Phjong, Ninh Binh Province	Viet Nam	10/1996	VNM01	21.181970	106.00161	VD Viet Nam V1	VD Viet Nam V1-8	cerana	1773	76	10.7	46623619	29994005	64.3	26399304	HiSeq 4000
V475_1	VD475_1	V. destructor	A. cerana	Hoc Chau	Viet Nam	11/1998	VNM02	20.922080	104.75209	VD Korea K1	VD Korea K1-1/K1-2	mellifera	6	555	19.3	50827244	48381793	95.2	44842843	HiSeq 4000
V153_2	VD153_2	V. destructor	A. mellifera	Jungup	South Korea	10/1996	KOR01	35.569880	126.85589	VD Korea K1	VD Korea K1-1/K1-2	mellifera	0	5824	22.2	56113114	55021332	98.1	51939126	HiSeq 4000
V159_1	VD159_1	V. destructor	A. mellifera	Jungup	South Korea	10/1996	KOR01	35.569880	126.85589	VD Korea K1	VD Korea K1-1/K1-2	mellifera	14	5707	24.2	64315412	60528716	94.1	55776173	HiSeq 4000
V639_1	VD639_1	V. destructor	A. mellifera	Heilongjiang Province	China	09/2001	CHN01	46.617160	131.15727	VD Korea K1	VD Korea K1-1/K1-2	mellifera	20	1096	17.2	45259621	43071591	95.2	39536470	HiSeq 4000
V625_2	VD625_2	V. destructor	A. mellifera	Fengman, Lishu Co., Changchun, Jilin Province	China	05/2001	CHN02	43.821600	126.56227	VD Korea K1	VD Korea K1-1/K1-2	mellifera	2	217	16.6	46030858	41845105	90.9	37822659	HiSeq 4000
V622_1	VD622_1	V. destructor	A. mellifera	Apiary of Jinzhong, Yuci, Shanxi Province	China	06/2001	CHN03	37.697790	112.70822	VD Korea K1	VD Korea K1-1/K1-2	mellifera	0	1243	13.8	39024603	35070629	89.9	31675460	HiSeq 4000
V641_1	VD641_1	V. destructor	A. mellifera	Mianzhu County, Sichuan Province	China	09/2001	CHN04	31.338070	104.22074	VD Korea K1	VD Korea K1-7	mellifera	6	871	21.0	54493356	52402973	96.2	48743391	HiSeq 4000
V646_1	VD646_1	V. destructor	A. mellifera	Zhongshan, Guangdong Province	China	11/2001	CHN09	23.927990	113.15883	VD Korea K1	VD Korea K1-7	mellifera	4	322	25.8	66768826	64203224	96.2	59914828	HiSeq 4000
V150_2	VD150_2	V. destructor	A. mellifera	Nho Quan	Viet Nam	10/1996	VNM03	20.322630	105.74960	VD Korea K1	VD Korea K1-1/K1-2	mellifera	26	218	33.2	86978010	84333771	97.0	80000889	HiSeq 4000
V474_1	VD474_1	V. destructor	A. mellifera	Hoc Chau	Viet Nam	01/1998	VNM02	20.922080	104.75209	VD Korea K1	VD Korea K1-1/K1-2	mellifera	52	648	24.8	63242091	61585346	97.4	58056817	HiSeq 4000
V028	VJ028	V. jacobsoni	A. cerana	UPM, Serdang, Kuala Lumpur	Malaysia	03/1995	MYS01	2.987890	101.73517	VJ Malaysia 2	VJ Malaysia 2-1	cerana	168	66	50.4	131769709	125151456	95.0	120488777	HiSeq 4000
V780	VJ780	V. jacobsoni	A. cerana	Chauthanh Village, Bentra District	Viet Nam	03/2003	VNM06	10.304830	106.35520	VJ Laos L1	VJ Laos L1-3	cerana	40372	4	15.6	47957642	38953995	81.2	36880992	HiSeq 4000
V787_2	VJ787_2	V. jacobsoni	A. cerana	Rach Gia Town, Kien Giang	Viet Nam	08/2004	VNM04	10.021500	105.09109	VJ Laos L1	VJ Laos L1-4	cerana	2325	28	29.6	78330554	75988942	97.0	72622136	HiSeq 4000
V788_1 (replicate discarded)	VJ788_1	V. jacobsoni	A. cerana	My tho City, Kien Giang	Viet Nam	08/2004	VNM05	10.376520	106.34388	VJ Laos L1	VJ Laos L1-4	cerana	1291	2	7.0	43611088	19916649	45.7	18505755	HiSeq 4000
V788_2	VJ788_2	V. jacobsoni	A. cerana	My tho City, Kien Giang	Viet Nam	08/2004	VNM05	10.376520	106.34388	VJ Laos L1	VJ Laos L1-4	cerana	72	2	10.8	30299193	27206968	89.8	25941038	HiSeq 4000
V789_1	VJ789_1	V. jacobsoni	A. cerana	Giong Trom Town, Ben Tre	Viet Nam	08/2004	VNM07	10.157180	106.50445	VJ Laos L1	VJ Laos L1-4	cerana	859	4	23.6	71558906	58473472	81.7	56649683	HiSeq 4000
V333_1	VJ333_1	V. jacobsoni	A. cerana	Parung Panjang, Java	Indonesia	06/1998	IDN01	-6.348140	106.56935	VJ Java 1	VJ Java 1-1	cerana	21	10	17.9	46019459	44634047	97.0	42592178	HiSeq 4000
V341_1	VJ341_1	V. jacobsoni	A. cerana	Sukabumi, West Java	Indonesia	06/1998	IDN02	-6.923870	106.93156	VJ Java 1	VJ Java 1-2	cerana	119	4	24.8	63057170	61695629	97.8	59251392	HiSeq 4000
V347_1	VJ347_1	V. jacobsoni	A. cerana	Cililin	Indonesia	06/1998	IDN04	-6.983450	107.43608	VJ Java 1	VJ Java 1-3	cerana	61	4	18.2	48579482	45809925	94.3	42898293	HiSeq 4000
V363_1	VJ363_1	V. jacobsoni	A. cerana	Johor, East Java	Indonesia	06/1998	IDN05	-7.223850	112.73320	VJ Java 1	VJ Java 1-4	cerana	5115	64	21.3	54026529	52825280	97.8	50723792	HiSeq 4000
V565_2	VJ565_2	V. jacobsoni	A. cerana	Gunung Arca, Sukabumi	Indonesia	04/2002	IDN03	-6.926560	106.95533	VJ Java 1	VJ Java 1-5	cerana	714	28	42.2	107959838	106090376	98.3	102687165	HiSeq 4000
V325_1	VJ325_1	V. jacobsoni	A. cerana	Kuta, Lombok	Indonesia	04/1998	IDN06	-8.880270	116.28361	VJ Lombok 2	VJ Lombok 2-1	cerana	152	0	24.7	62599084	61346630	98.0	58922351	HiSeq 4000
V950_1	VJ950_1	V. jacobsoni	A. cerana	Aiyura	Papua New Guinea	05/1997	PNG01	-6.343340	145.90809	VJ PNG 1	VJ PNG 1-1	cerana	475	68	25.8	66335209	64182154	96.8	61215275	HiSeq 4000
V854_1	VJ854_1	V. jacobsoni	A. cerana	Goroka	Papua New Guinea	05/2008	PNG04	-6.083470	145.38626	VJ PNG 1	VJ PNG 1-1	cerana	5230	12	16.7	44613684	42039605	94.2	39498649	HiSeq 4000
V853_4	VJ853_4	V. jacobsoni	A. cerana	Daru	Papua New Guinea	06/2008	PNG11	-9.078190	143.20997	VJ PNG 1	VJ PNG 1-2	cerana	200	60	10.4	35177854	32209884	91.6	29491225	HiSeq 4000
V983_1	VJ983_1	V. jacobsoni	A. cerana	Asaro	Papua New Guinea	11/2015	PNG05	-6.012000	145.32500	VJ PNG 1	VJ PNG 1-1	cerana	284	102	19.1	50283170	47881277	95.2	44989750	HiSeq 4000
V847	VJ847	V. jacobsoni	A. mellifera	Wabag, Enga Province	Papua New Guinea	06/2008	PNG10	-5.482670	143.66373	VJ PNG 1	VJ PNG 1-1	mellifera	23	115346	22.5	64504980	56218551	87.2	52755225	HiSeq 4000
V852_3	VJ852_3	V. jacobsoni	A. mellifera	Goroka, Eastern Highlands Province	Papua New Guinea	05/2008	PNG04	-6.083470	145.38626	VJ PNG 1	VJ PNG 1-1	mellifera	4	2288	20.9	55816067	54523224	97.7	51390551	HiSeq 4000
V857_2	VJ857_2	V. jacobsoni	A. mellifera	Benabena, Eastern Highlands Province	Papua New Guinea	05/2008	PNG03	-6.128125	145.42904	VJ PNG 1	VJ PNG 1-1	mellifera	20	429	18.8	57477057	46759609	81.4	44551404	HiSeq 4000
V856_1	VJ856_1	V. jacobsoni	A. mellifera	Henganofi Border, Eastern Highlands Province	Papua New Guinea	05/2008	PNG02	-6.260350	145.59681	VJ PNG 1	VJ PNG 1-1	mellifera	15	9803	13.6	58525652	37045620	63.3	31803798	HiSeq 4000
V952_1	VJ952_1	V. jacobsoni	A. mellifera	Amayufa, Eastern Highlands Province	Papua New Guinea	05/2014	PNG06	-5.975581	145.27676	VJ PNG 1	VJ PNG 1-1	mellifera	18	36	22.3	61969612	55519034	89.6	52870625	HiSeq 4000
V953_2	VJ953_2	V. jacobsoni	A. mellifera	Benabena, Eastern Highlands Province	Papua New Guinea	05/2014	PNG03	-6.128125	145.42904	VJ PNG 1	VJ PNG 1-1	mellifera	22	26196	22.3	62054570	56103596	90.4	53032044	HiSeq 4000
V954_2	VJ954_2	V. jacobsoni	A. mellifera	Windy Lodge	Papua New Guinea	06/2014	PNG07	-6.035146	144.96032	VJ PNG 1	VJ PNG 1-1	mellifera	14	159	22.5	58553267	56180797	96.0	52846232	HiSeq 4000
V955_1	VJ955_1	V. jacobsoni	A. mellifera	Donakanage	Papua New Guinea	06/2014	PNG08	-5.882240	144.82077	VJ PNG 1	VJ PNG 1-1	mellifera	15	460	24.7	64865555	61577522	94.9	58553870	HiSeq 4000
V956_4	VJ956_4	V. jacobsoni	A. mellifera	Minj	Papua New Guinea	06/2014	PNG09	-5.826640	144.40741	VJ PNG 1	VJ PNG 1-2	cerana	1444	69	17.1	47824500	45966716	96.1	42702569	HiSeq 4000
Mom #1	Mom #1	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	02/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	0	7448	20.2	52944228	50821794	96.0	52944228	HiSeq 4000
Son #1	Son #1	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	02/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	0	341	19.8	51620028	49266370	95.4	51620028	HiSeq 4000
Mom #2	Mom #2	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	02/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	0	50765	24.8	65220912	62435130	95.7	65220912	HiSeq 4000
Son #2	Son #2	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	02/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	0	422	18.7	48499085	46511492	95.9	48499085	HiSeq 4000
Mom #7	Mom #7	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	02/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	0	563	16.3	43293860	41401282	95.6	43293860	HiSeq 4000
Son #7	Son #7	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	02/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	0	660	18.8	49134877	46899128	95.5	49134877	HiSeq 4000
Mom #15	Mom #15	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	05/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	0	100	22.5	59835146	56138955	93.8	59835146	HiSeq 4000
Son #15	Son #15	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	05/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	1	193	21.4	57025270	53922121	94.6	57025270	HiSeq 4000
Mom #17	Mom #17	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	05/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	2	40283	22.1	58014781	55352009	95.4	58014781	HiSeq 4000
Son #17	Son #17	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	05/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	2	949	33.0	87139528	82953771	95.2	87139528	HiSeq 4000
Mom #19	Mom #19	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	05/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	0	24181	18.1	47748745	45343415	95.0	47748745	HiSeq 4000
Son #19	Son #19	V. destructor	A. mellifera	OIST Ecology and Evolution experimental apiary, Okinawa	Japan	05/2018	JPN01	26.456448	127.82386	VD Korea K1	VD Korea K1-1/K1-2	mellifera	5	498	24.5	65304925	61543197	94.2	65304925	HiSeq 4000

  #row_spec(1:23, bold = T, color = "white", background = "#E60000") %>%
  #row_spec(24:39, bold = T, color = "white", background = "#FF9999")  %>%
  #row_spec(40:74, bold = T, color = "white", background = "#0000FF")  %>%
  #row_spec(75:96, bold = T, color = "white", background = "#3399FF")

# Download the country borders layer
#download.file("http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip" , destfile="world_shape_file.zip")
#system("unzip world_shape_file.zip")
world_spdf=readOGR(dsn= getwd() , layer="TM_WORLD_BORDERS_SIMPL-0.3")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/alphamanae/Documents/GitHub/varroa-host-jump/R_data", layer: "TM_WORLD_BORDERS_SIMPL-0.3"
## with 246 features
## It has 11 fields
## Integer64 fields read as strings:  POP2005

data.vdac <- metadata %>% filter(Species == "V. destructor") %>% filter(Host_collected == "A. cerana") 

data.vdam <- metadata %>% filter(Species == "V. destructor") %>% filter(Host_collected == "A. mellifera")

data.vjac <- metadata %>% filter(Species == "V. jacobsoni") %>% filter(Host_collected == "A. cerana") 

data.vjam <- metadata %>% filter(Species == "V. jacobsoni") %>% filter(Host_collected == "A. mellifera") 

The following interactive map showed the distribution of both V. destructor and V. jacobsoni samples on either their original host the eastern honey bee A. cerana or on the new host the western honey bee A. mellifera.

If you pass your mouse over a particular point, a pop-up will appear presenting essential information about the particular Varroa mite sample. On the right top corner, a layer control button allow to subselect species/host couple.

V. destructor on A. cerana
V. destructor on A. mellifera
V. jacobsoni on A. cerana
V. jacobsoni on A. mellifera

## Create the interactive map with leaflet
leaflet(metadata) %>% 
  addTiles(group = "OSM (default)") %>%
  ## add the layer other than default we would like to use for background
  addProviderTiles(providers$CartoDB.PositronNoLabels, group = "Positron NoLabels") %>%
  ## adding the layers with V. destructor
  addCircleMarkers(data = data.vdac, data.vdac$Approximate.Y, data.vdac$Approximate.X,
                   weight = 0.5,
                   col = "#FB0000", 
                   radius = 4, 
                   fillOpacity = 0.9, 
                   stroke = T, 
                   popup = ~paste(sep = "<br/>",
                              paste ("Name: ", as.character(ID)),
                              paste ("Host: ", as.character(Host_collected)),
                              paste ("Date: ", as.character(Date)),
                              paste ("mtDNA: ", as.character(Haplotype))),
                   group = 'V. destructor on A. cerana') %>%

  addCircleMarkers(data = data.vdam, data.vdam$Approximate.Y, data.vdam$Approximate.X,
                   weight = 0.5,
                   col = "#FF9999", 
                   radius = 4, 
                   fillOpacity = 0.9, 
                   stroke = T, 
                   popup = ~paste(sep = "<br/>",
                              paste ("Name: ", as.character(ID)),
                              paste ("Host: ", as.character(Host_collected)),
                              paste ("Date: ", as.character(Date)),
                              paste ("mtDNA: ", as.character(Haplotype))),
                   group = 'V. destructor on A. mellifera') %>%

  ## adding the layers with V. jacobsoni
  addCircleMarkers(data = data.vjac, data.vjac$Approximate.Y, data.vjac$Approximate.X,
                   weight = 0.5,
                   col = "#0000FF", 
                   radius = 4, 
                   fillOpacity = 0.9, 
                   stroke = T, 
                   popup = ~paste(sep = "<br/>",
                              paste ("Name: ", as.character(ID)),
                              paste ("Host: ", as.character(Host_collected)),
                              paste ("Date: ", as.character(Date)),
                              paste ("mtDNA: ", as.character(Haplotype))),
                   group = 'V. jacobsoni on A. cerana') %>%

  addCircleMarkers(data = data.vjam, data.vjam$Approximate.Y, data.vjam$Approximate.X,
                   weight = 0.5,
                   col = "#00CCFF", 
                   radius = 4, 
                   fillOpacity = 0.9, 
                   stroke = T,    
                   popup = ~paste(sep = "<br/>",
                              paste ("Name: ", as.character(ID)),
                              paste ("Host: ", as.character(Host_collected)),
                              paste ("Date: ", as.character(Date)),
                              paste ("mtDNA: ", as.character(Haplotype))),
                   group = 'V. jacobsoni on A. mellifera') %>%

    ## adding the control button to remove or add layers of points 
  addLayersControl(position = "topright",
    baseGroups = c("OSM (default)", "Positron NoLabels"),
    overlayGroups = c("V. destructor on A. cerana", 
                      "V. destructor on A. mellifera", 
                      "V. jacobsoni on A. cerana", 
                      "V. jacobsoni on A. mellifera"), 
    options = layersControlOptions(collapsed = TRUE))  %>%
  ## show the positron background prerably to the OSM layer
  showGroup("Positron NoLabels")

2. Reads and mapping statistics for Varroa libraries

Here we plot the mean read depth computed from each .bam file mapped to the reference vdes3.0.fasta using samtools depth -a FILE.bam | awk '{c++;s+=$3}END{print s/c}' in relation to the mapping percentage.

As you can see by moving your pointer on the graph, the individual V788_1 has been removed from the analysis and replaced by the replicate V788_2 with much better reads statistics.

raichu <- ggplot(metadata, aes(x=Mapping_rate, y = Mean_Read_Depth, col=Species, shape=Host_collected, label = ID))
mycol = c("red2", "blue2")
raichu <- raichu + geom_point(size = 3)
raichu <- raichu + scale_color_manual(values = mycol)
raichu <- raichu + scale_size_manual(values = mycol)
raichu <- raichu +  theme_calc()
raichu <- raichu +  xlab("Mapping rate to Vdes_3.0 reference genome(%)")
raichu <- raichu +  ylab("Mean average read depth")
raichu <- raichu +  ylim(0,20)
#raichu

## make an interactive version of the scatter plot
intraichu <- ggplotly(raichu)
intraichu

We highlight the mean read depth regarding the general distribution, where is each group of individuals “species per host” located in it.

evoli <- ggplot(metadata, aes(Mean_Read_Depth, fill = Species))
evoli <- evoli + geom_histogram(bins = 44) 
evoli <- evoli + gghighlight() 
evoli <- evoli + facet_wrap(~Species~Host_collected)
evoli <- evoli + scale_fill_manual(values = c("red", "blue"))
evoli <- evoli + theme_classic()
evoli <- evoli + labs(x = "Mean Read Depth using NextGenMap", y = "Number of samples")
evoli <- evoli + ggtitle("Read depth does not appear as unbalanced accross species and host")
evoli

Following is a plot of the total number of reads generated per individual and the proportion of reads mapped onto the reference.

togepi <- ggplot(metadata, aes(x=ID, label = Mapping_rate))
togepi <- togepi + geom_bar(aes(weight=Total_reads), fill="black", position="dodge")
togepi <- togepi + geom_bar(aes(weight=Mapped), fill="#1d96f2", position="dodge")
togepi <- togepi + theme(axis.text.x = element_text(angle = 90, hjust = 1))
togepi <- togepi + labs(x = "Individuals", y = "Number of reads")
togepi <- togepi + ggtitle("Proportion of reads mapped against reference using NextGenMap")
togepi

3. Variation in the body size

Before proceeding to DNA extraction with a destructive method, we measured the body size of each Varroa mite whenever possible on both dorsal and ventral view. We measured seven morphological characters but retained only the body width and length for comparison here between Varroa species and their associated honey bee hosts.

The scatterplot also shows the average values obtained by Anderson and Trueman (2000) Exp. App. Acar. for V. destructor, V. jacobsoni and undetermined mites from Luzon and Mindanao.

# Table with 512 morpho data for the moment
morpho <- read.csv("Varroa-morpho-subset.csv", header = TRUE)
#kable(cbind(morpho)) %>%
#  kable_styling(bootstrap_options = "striped", font_size = 10) %>%
#  scroll_box(width = "100%", height = "400px")

morpho %>% 
  plot_ly(x = ~BW_V, y = ~BL_V,
          color = ~Species,
          symbol = ~Host,
          hoverinfo = "text",
          text = ~paste("Name: ", Sample.ID, "<br>",
                        "Species: ", Species, "<br>",
                        "Host: ", Host, "<br>",
                        "Country: ", Country.Island)) %>%
  add_markers(colors = c("indianred2", "red3", "skyblue2", "blue2", "black"), 
              symbols = c( "triangle-up", "circle", "x"), 
              size = 3) %>%
  layout(xaxis = list(title = "Body width (mm)"),
         yaxis = list(title = "Body length (mm)"),
         title = "Variation in body size measured from ventral view")

morpho %>% 
  plot_ly(x = ~BW_D, y = ~BL_D,
          color = ~Species,
          symbol = ~Host,
          hoverinfo = "text",
          text = ~paste("Name: ", Sample.ID, "<br>",
                        "Species: ", Species, "<br>",
                        "Host: ", Host, "<br>",
                        "Country: ", Country.Island)) %>%
  add_markers(colors = c("indianred2", "red3", "skyblue2", "blue2", "black"), 
              symbols = c( "triangle-up", "circle", "x"), 
              size = 3) %>%
  layout(xaxis = list(title = "Body width (mm)"),
         yaxis = list(title = "Body length (mm)"),
         title = "Variation in body size measured from dorsal view")

STATS TEST?

4. Genetic diversity on 43 samples and two species confonded

2D and 3D PCAs on ~ 1.4 million filtered biallelic SNPs

After formatting the variant call file .vcf into a plink format by renaming each chromosome (e.g., NW_019211454.1 by 1, NW_019211455.1 by 2, …), and applying the option pca, we obtained the two files .eigenval and .eigenvec.

The following 2D and 3D PCAs are based on the ~1.4 million biallelic SNPs with a minor allelic frequency higher than 5%.

pca.all<- read_table2("43indplink.eigenvec", col_names = FALSE)
eigenval.all <- scan("43indplink.eigenval")

## Here the first two columns contains the individuals names as we did not specify it in plink earlier
# We remove the first column and assign the new name of column 1
pca.all <- pca.all[,-1]
names(pca.all)[1] <- "ID_NAMES"

# Rename the columns for PCA axis 
names(pca.all)[2:ncol(pca.all)] <- paste0("PC", 1:(ncol(pca.all)-1))

## Joint both table by using the FILE_ID as common parameters
metapca.all <- left_join(pca.all, metadata, by = c("ID_NAMES" = "SEQ_ID"))

### PLOTS
# first convert to percentage variance explained
pve <- data.frame(PC = 1:20, pve = eigenval.all/sum(eigenval.all)*100)

# make Eigenplot
a <- ggplot(pve, aes(PC, pve)) + geom_bar(stat = "identity")
a + ylab("Percentage variance explained") + theme_light()

metapca.all %>% 
  plot_ly(x = ~PC1, y = ~PC2,
          color = ~Species,
          symbol = ~Host_collected,
          hoverinfo = "text",
          text = ~paste("Name: ", ID, "<br>",
                        "Species: ", Species, "<br>",
                        "Host: ", Host_collected, "<br>",
                        "Country: ", Country)) %>%
  add_markers(colors = c("red3", "blue2"), 
              symbols = c( "triangle-up", "circle"), 
              size = 4)%>%
  layout(title = "Varroa mites genetically differentiate between hosts in both species")

metapca.all %>% 
  plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
          color = ~Country,
          hoverinfo = "text",
          text = ~paste("Name: ", ID, "<br>",
                        "Species: ", Species, "<br>",
                        "Host: ", Host_collected, "<br>",
                        "Country: ", Country,
                        "Haplogroup: ", COI_haplogroup)) %>%
  add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"), 
              size = 10) %>%
  layout(title = "Mites on A. cerana genetically differs within the native range")

metapca.all %>% 
  plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
          color = ~COI_haplogroup,
          hoverinfo = "text",
          text = ~paste("Name: ", ID, "<br>",
                        "Species: ", Species, "<br>",
                        "Host: ", Host_collected, "<br>",
                        "Country: ", Country,
                        "Haplogroup: ", COI_haplogroup)) %>%
  add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"), 
              size = 10) %>%
  layout(title = "Mitochondrial lineages can be a pre-indicator of whole genome differentiation")

2D and 3D PCAs on ~ 200,000 biallelic SNPs pruned for LD

Same idea as the previous PCA, except that only 236,713 SNPs were kept to reduce the effect of linkage disequilibrium.

We can see that even with a smaller reduce number of SNPs, Varroa mites species are genetically distinct and

pca.allLD<- read_table2("43ind-ldpruned.eigenvec", col_names = FALSE)
eigenval.allLD <- scan("43ind-ldpruned.eigenval")

## Here the first two columns contains the individuals names as we did not specify it in plink earlier
# We remove the first column and assign the new name of column 1
pca.allLD <- pca.allLD[,-1]
names(pca.allLD)[1] <- "ID_NAMES"

# Rename the columns for PCA axis 
names(pca.allLD)[2:ncol(pca.allLD)] <- paste0("PC", 1:(ncol(pca.allLD)-1))

## Joint both table by using the FILE_ID as common parameters
metapca.allLD <- left_join(pca.allLD, metadata, by = c("ID_NAMES" = "SEQ_ID"))

### PLOTS
# first convert to percentage variance explained
pve <- data.frame(PC = 1:20, pve = eigenval.allLD/sum(eigenval.allLD)*100)

# make Eigenplot
a <- ggplot(pve, aes(PC, pve)) + geom_bar(stat = "identity")
a + ylab("Percentage variance explained") + theme_light()

metapca.allLD %>% 
  plot_ly(x = ~PC1, y = ~PC2,
          color = ~Species,
          symbol = ~Host_collected,
          hoverinfo = "text",
          text = ~paste("Name: ", ID, "<br>",
                        "Species: ", Species, "<br>",
                        "Host: ", Host_collected, "<br>",
                        "Country: ", Country)) %>%
  add_markers(colors = c("red3", "blue2"), 
              symbols = c( "triangle-up", "circle"), 
              size = 4) %>%
  layout(title = "~200,000 LD pruned SNPs still allowed to see host genetic divergence")

metapca.allLD %>% 
  plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
          color = ~Country,
          hoverinfo = "text",
          text = ~paste("Name: ", ID, "<br>",
                        "Species: ", Species, "<br>",
                        "Host: ", Host_collected, "<br>",
                        "Country: ", Country,
                        "Haplogroup: ", COI_haplogroup)) %>%
  add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"), 
              size = 10) %>%
  layout(title = "3D PCA color-coded by country")

metapca.allLD %>% 
  plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
          color = ~COI_haplogroup,
          hoverinfo = "text",
          text = ~paste("Name: ", ID, "<br>",
                        "Species: ", Species, "<br>",
                        "Host: ", Host_collected, "<br>",
                        "Country: ", Country,
                        "Haplogroup: ", COI_haplogroup)) %>%
  add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"), 
              size = 10) %>%
  layout(title = "3D PCA color-coded by mtDNA COX1 haplogroup")

5. Demography and best scenario bootstraps

After estimating the demographic history for each Varroa species on their novel and original hosts using fastsimcoal2, we plotted here the distribution of the log likelihood from the 100 run replicates for each scenario and each SFS subset size.

This method for model selection comes on top of the AIC calculation and is followed the same methods as described in Meier et al. 2016. Mol. Ecol.. For the AIC values, see Table S7.

Varroa destructor demographic history

Likelihood distribution for the four VD SFS datasets

vd.lhoods1050 <- read.csv("fsc_run/lhoods_1050SNPS_vdno475.csv", header = TRUE)
vd.lhoods10500 <- read.csv("fsc_run/lhoods_10500SNPS_vdno475.csv", header = TRUE)
vd.lhoods35000 <- read.csv("fsc_run/lhoods_35000SNPS_vdno475.csv", header = TRUE)
vd.lhoods105000 <- read.csv("fsc_run/lhoods_105000SNPS_vdno475.csv", header = TRUE)

par(mfrow=c(1,4))
#Plot 1 for likelihoods with smaller subset
boxplot(range = 0,vd.lhoods1050$scenario1_div [3:102], vd.lhoods1050$scenario2_divGwt2n[3:102], vd.lhoods1050$scenario3_divGwt[3:102], vd.lhoods1050$scenario4_IM[3:102], vd.lhoods1050$scenario5_IMGwt2n[3:102], vd.lhoods1050$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vd.lhoods1050$scenario1_div [2],col = "red", pch = 20)
points(2,vd.lhoods1050$scenario2_divGwt2n [2],col = "red", pch = 20)
points(3,vd.lhoods1050$scenario3_divGwt [2],col = "red", pch = 20)
points(4,vd.lhoods1050$scenario4_IM [2],col = "red", pch = 20)
points(5,vd.lhoods1050$scenario5_IMGwt2n [2],col = "red", pch = 20)
points(6,vd.lhoods1050$scenario6_IMGwt [2],col = "red", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. destructor 377 SNPS")

#Plot 2 for likelihoods with second smaller subset
boxplot(range = 0,vd.lhoods10500$scenario1_div [3:102], vd.lhoods10500$scenario2_divGwt2n[3:102], vd.lhoods10500$scenario3_divGwt[3:102], vd.lhoods10500$scenario4_IM[3:102], vd.lhoods10500$scenario5_IMGwt2n[3:102], vd.lhoods10500$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vd.lhoods10500$scenario1_div [2],col = "red", pch = 20)
points(2,vd.lhoods10500$scenario2_divGwt2n [2],col = "red", pch = 20)
points(3,vd.lhoods10500$scenario3_divGwt [2],col = "red", pch = 20)
points(4,vd.lhoods10500$scenario4_IM [2],col = "red", pch = 20)
points(5,vd.lhoods10500$scenario5_IMGwt2n [2],col = "red", pch = 20)
points(6,vd.lhoods10500$scenario6_IMGwt [2],col = "red", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. destructor 3,760 SNPS")

#Plot 3 for likelihoods with third subset
boxplot(range = 0,vd.lhoods35000$scenario1_div [3:102], vd.lhoods35000$scenario2_divGwt2n[3:102], vd.lhoods35000$scenario3_divGwt[3:102], vd.lhoods35000$scenario4_IM[3:102], vd.lhoods35000$scenario5_IMGwt2n[3:102], vd.lhoods35000$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vd.lhoods35000$scenario1_div [2],col = "red", pch = 20)
points(2,vd.lhoods35000$scenario2_divGwt2n [2],col = "red", pch = 20)
points(3,vd.lhoods35000$scenario3_divGwt [2],col = "red", pch = 20)
points(4,vd.lhoods35000$scenario4_IM [2],col = "red", pch = 20)
points(5,vd.lhoods35000$scenario5_IMGwt2n [2],col = "red", pch = 20)
points(6,vd.lhoods35000$scenario6_IMGwt [2],col = "red", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. destructor 12,665 SNPS")

#Plot 4 for likelihoods with larger subset
boxplot(range = 0,vd.lhoods105000$scenario1_div [3:102], vd.lhoods105000$scenario2_divGwt2n[3:102], vd.lhoods105000$scenario3_divGwt[3:102], vd.lhoods105000$scenario4_IM[3:102], vd.lhoods105000$scenario5_IMGwt2n[3:102], vd.lhoods105000$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vd.lhoods105000$scenario1_div [2],col = "red", pch = 20)
points(2,vd.lhoods105000$scenario2_divGwt2n [2],col = "red", pch = 20)
points(3,vd.lhoods105000$scenario3_divGwt [2],col = "red", pch = 20)
points(4,vd.lhoods105000$scenario4_IM [2],col = "red", pch = 20)
points(5,vd.lhoods105000$scenario5_IMGwt2n [2],col = "red", pch = 20)
points(6,vd.lhoods105000$scenario6_IMGwt [2],col = "red", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. destructor 38,108 SNPS")

After summarizing the 100 best replicates from the 100 simulated SFS run, we computed the mean estimate value and 95% confidence interval for each demographic parameter using boot package.

Scenario 6: IM + BOT bootstraps computation

my.mean = function(x, indices) {
return( mean( x[indices] ) )
}

# import the file containing the 100 best runs
vd.boot6 <- read.csv("fsc_run/vdno475_105000nopv_BOOT.csv", header = TRUE)

# Effective size for Varroa mites on original host Apis cerana
N_VdAc1.boot = boot(vd.boot6$NVDAC1, my.mean, 10000)
my.mean(vd.boot6$NVDAC1,1:length(vd.boot6$NVDAC1))

## [1] 412.5

boot.ci(N_VdAc1.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = N_VdAc1.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (400.6, 424.4 )   (400.7, 424.1 )  
## 
## Level     Percentile            BCa          
## 95%   (400.9, 424.3 )   (400.8, 424.3 )  
## Calculations and Intervals on Original Scale

# Effective size for Varroa mites on new host Apis mellifera
N_VdAm0.boot = boot(vd.boot6$NVDAM0, my.mean, 10000)
my.mean(vd.boot6$NVDAM0,1:length(vd.boot6$NVDAM0))

## [1] 568.11

boot.ci(N_VdAm0.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = N_VdAm0.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (550.2, 586.2 )   (550.4, 586.6 )  
## 
## Level     Percentile            BCa          
## 95%   (549.6, 585.8 )   (549.5, 585.7 )  
## Calculations and Intervals on Original Scale

# Number of generations when new population of Varroa mites was founded
TBOT.boot = boot(vd.boot6$TBOT, my.mean, 10000)
my.mean(vd.boot6$TBOT,1:length(vd.boot6$TBOT))

## [1] 722.98

boot.ci(TBOT.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = TBOT.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (690.0, 755.8 )   (689.9, 755.7 )  
## 
## Level     Percentile            BCa          
## 95%   (690.3, 756.1 )   (690.4, 756.1 )  
## Calculations and Intervals on Original Scale

# Growth rate of Varroa mites on new host Apis mellifera
GAM.boot = boot(vd.boot6$GAM, my.mean, 10000)
my.mean(vd.boot6$GAM,1:length(vd.boot6$GAM))

## [1] -0.0034098

boot.ci(GAM.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = GAM.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (-0.0036, -0.0032 )   (-0.0036, -0.0032 )  
## 
## Level     Percentile            BCa          
## 95%   (-0.0036, -0.0032 )   (-0.0036, -0.0032 )  
## Calculations and Intervals on Original Scale

# Effective size for Varroa mites able to jump on A. mellifera to create new population
JUMP.boot = boot(vd.boot6$JUMP, my.mean, 10000)
my.mean(vd.boot6$JUMP,1:length(vd.boot6$JUMP))

## [1] 165.04

boot.ci(JUMP.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = JUMP.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (160.3, 169.7 )   (160.3, 169.6 )  
## 
## Level     Percentile            BCa          
## 95%   (160.5, 169.8 )   (160.5, 169.8 )  
## Calculations and Intervals on Original Scale

N0M1.boot = boot(vd.boot6$N0M1, my.mean, 10000)
my.mean(vd.boot6$N0M1,1:length(vd.boot6$N0M1))

## [1] 0.0005941365

boot.ci(N0M1.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = N0M1.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0.0003,  0.0009 )   ( 0.0003,  0.0008 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0.0003,  0.0009 )   ( 0.0004,  0.0011 )  
## Calculations and Intervals on Original Scale

N1M0.boot = boot(vd.boot6$N1M0, my.mean, 10000)
my.mean(vd.boot6$N1M0,1:length(vd.boot6$N1M0))

## [1] 0.7604321

boot.ci(N1M0.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = N1M0.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0.7586,  0.7623 )   ( 0.7586,  0.7623 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0.7586,  0.7623 )   ( 0.7586,  0.7623 )  
## Calculations and Intervals on Original Scale

MACtoAM.boot = boot(vd.boot6$MACtoAM, my.mean, 10000)
my.mean(vd.boot6$MACtoAM,1:length(vd.boot6$MACtoAM))

## [1] 0.001885131

boot.ci(MACtoAM.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = MACtoAM.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0.0018,  0.0019 )   ( 0.0018,  0.0019 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0.0018,  0.0019 )   ( 0.0018,  0.0019 )  
## Calculations and Intervals on Original Scale

MAMtoAC.boot = boot(vd.boot6$MAMtoAC, my.mean, 10000)
my.mean(vd.boot6$MAMtoAC,1:length(vd.boot6$MAMtoAC))

## [1] 1.044978e-06

boot.ci(MAMtoAC.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = MAMtoAC.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0,  0 )   ( 0,  0 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0,  0 )   ( 0,  0 )  
## Calculations and Intervals on Original Scale

Here we modified the .par file from the best scenario 6 for V. destructor by adding the mean values computed in the step before. Using the ParFileInterpreter-v6.3.1.r developed by Victor Sousa and available on fastsimcoal2 website, we can now visulized the best demographic scenario.

Draw parameters for best scenario with mean bootstraps

#...........................................................................................
# (c) Laurent Excoffier and Vitor Sousa June-November 2015
# 
# Small R program to draw the evolutonary scenario described by a given par file
# This is mainly for visual checking that the modeled scenario corresponds to   
# what was intended
#
#...........................................................................................
#
args=commandArgs(TRUE)
print(args)

## character(0)

#Expects par file name on command line

if(length(args)) {
  parFileName=args[1]
} else {  
  #REPLACE BY THE NAME OF THE PAR FILE YOU WANT TO ANALYSE
  parFileName="fsc_run/VDM6_downsample_IMGwt_49_boot.par"
}

separator=" "
migrMatCol="coral"
admixCol="blue"
popFusionColor="black"
popCol="lightgrey"
popBorderCol="black"
inbreedColor="lightblue"
oldSampColor="olivedrab4"
timeCol="tan4"
growthCol="hotpink3"
propLastsegment=0.05
migMatNameProp=0.8
migMatLineLength=0.3
timeProp=0.6
maxRadius=1/40
minRadius=maxRadius/3
arrowLength=0.2
logScaleAxis=""
timeOffset=0.25
migrOffset=0.1
curvedArrowLTY=1
drawLogPopSize=T
plotMigrRates=T
migrRateTextSizeCEX=0.5

#Define plot area size for PDF
pdf.x.size=7
pdf.y.size=10

rescalingFactor=1.0 #Don't touch this!


printPDF=F

# parFile=readLines(con=parFileName)

###########################      Reading par file      #########################

parFile=scan(parFileName, character(0), sep = "\n", strip.white = TRUE) # separate each line

#pdfFileName=paste(parFileName, ".pdf", sep="")

#--- Clean input file by removing consecutive separators, and keep ---------------

#--- Function to remove separators within a string
removeTrailingSep=function(string, sep) {
  temp=strsplit(string, split=sep)
  temp2=temp[[1]][nchar(temp[[1]])>0]
  cleanStr=temp2[1]
  if (length(temp2)>1) {
    for (i in 2:length(temp2)) {
      cleanStr=paste(cleanStr, temp2[i], sep=sep)
    }
  }
  cleanStr
}

#--- Replace Keep by -9999
replaceKeep=function(string) {
  if (grepl("keep", string)) {
    return(gsub("keep", "-9999", string))
  }
  return(string)
}

#Remove both multiple consecutive whitespace and tabs and replace keep keyword
for (i in 1:length(parFile)) {
  parFile[i]=removeTrailingSep(parFile[i], sep='\t')
  parFile[i]=removeTrailingSep(parFile[i], sep=' ')
  parFile[i]=replaceKeep(parFile[i])
}

#-------------------------------------------------------------------------------

#--- Get number of samples on line 2 -----
l.numsamples=parFile[2]
sp.l=unlist(strsplit(l.numsamples, split=' '))
tab.l=unlist(strsplit(l.numsamples, split='\t'))
if (length(sp.l)>=length(tab.l)) {
  numSamples=as.numeric(sp.l[1])
  separator=" "
} else {
  numSamples=as.numeric(tab.l[1])
  separator="\t"
}

#--- Reading numbers on separate lines -----
getNumbers=function(start, parFile, numSamples) {  
  for (i in 1:numSamples) {
    curnum=as.numeric(unlist(strsplit(parFile[start+i], split=separator))[1])
    if (i==1) {
      num=curnum 
    } else {
      num=c(num, curnum)
    }
  }
  num
}

#--- Get population sizes -----------------
start=3
popSizes=getNumbers(start, parFile, numSamples)
#Rescaling pop sizes
popSizes=round(popSizes*rescalingFactor, digits=0)

iniPopSizes=popSizes

#--- Get sample sizes -----------------

readSampleSizesTimesAndInbreedingLevel=function(start, parFile, numSamples) {
  for (i in 1:numSamples) {
    curLine=unlist(strsplit(parFile[start+i], split=separator))
    curSampSize=as.numeric(curLine[1])
    curSampTime=0
    curInbreeding=0
    if (length(curLine)>1) curSampTime=as.numeric(curLine[2])
    if (length(curLine)>2) curInbreeding=as.numeric(curLine[3])
    if (i==1) {
      sampSize=curSampSize
      sampTime=curSampTime
      inbreeding=curInbreeding
    } else {
      sampSize=c(sampSize,curSampSize)
      sampTime=c(sampTime,curSampTime)
      inbreeding=c(inbreeding,curInbreeding)
    }
  }
  list(ss=sampSize, st=sampTime, inb=inbreeding)
}

start=start+numSamples+1
# sampSizes=getNumbers(start, parFile, numSamples)
sampSizesStats=readSampleSizesTimesAndInbreedingLevel(start, parFile, numSamples)
#Rescaling sample times
sampSizesStats$st=round(sampSizesStats$st*rescalingFactor, digits=0)

sampSizes=sampSizesStats$ss
sampTimes=sampSizesStats$st
inbrCoeff=sampSizesStats$inb

#--- Get growth rates -----------------
start=start+numSamples+1
growthRatesInitial=getNumbers(start, parFile, numSamples)
# save this into growthRates which will be used and updated when printing historical events
growthRates=growthRatesInitial

#--- Get number of migration matrices -----------------
start=start+numSamples+1
numMigMat=as.numeric(unlist(strsplit(parFile[start+1], split=separator))[1])

#--- Read migration matrix
readMigMat=function(start, parFile, numSamples) {
  for (i in 1:numSamples) {
    curmigs=as.numeric(unlist(strsplit(parFile[start+i], split=separator)))
    if (i==1) {
      migs=curmigs 
    } else {
      migs=rbind(migs, curmigs)
    }
  }
  rownames(migs)=1:numSamples
  migs 
}

#--- Get migration matrices as a list --------------
start=start+2
migMats=list()
if (numMigMat>0) {
  for (i in 1:numMigMat) {  
    curMigMat=readMigMat(start, parFile, numSamples) 
    migMats[[i]]=curMigMat
    start=start+numSamples+1
  }
}

#Rescaling migration rates
if (numMigMat>0) {
  for (i in 1:numMigMat) {
    migMats[[i]]=migMats[[i]]/rescalingFactor
  }
}

#--- Get number of historical events
start=start+1
numHistEvents=as.numeric(unlist(strsplit(parFile[start], split=separator))[1])

###################### HISTORICAL EVENTS HANDLING ##############################

#..... Read Historical Event .......
last.he.time=0
if (numHistEvents>0) {
  for (i in 1:numHistEvents) {
    start=start+1
    #Take care of nomig keyword
    nomig=F
    if (grepl("nomig", parFile[start])) {
      nomig=T
      gsub("nomig", "", parFile[start])
    }
    curHE=as.numeric(unlist(strsplit(parFile[start], split=separator)))
    if (nomig) curHE[7]=-1 
    #Rescaling time of event
    curHE[1]=round(curHE[1]*rescalingFactor, digits=0)
    if (i==1) {
      histEvents=curHE
      last.he.time=curHE[1]
    } else {
      histEvents=rbind(histEvents, curHE)
      if (histEvents[i,1] > last.he.time) last.he.time=histEvents[i,1] 
    }
  }
  if (numHistEvents>1) rownames(histEvents)=1:numHistEvents
}
names(last.he.time)=""
last.he.time=as.numeric(last.he.time)

yTimeLimit=0
if (last.he.time!=0) yTimeLimit=last.he.time*(1+propLastsegment) #Add propLastsegment to y axis after last event (to draw stuff)

#Reorder events by their times
if (numHistEvents>1) histEvents=histEvents[order(histEvents[,1],decreasing=FALSE),] else {
  if (numHistEvents==1) histEvents=matrix(histEvents, nrow=1,  byrow = T)
}

endReadParFile=start



##############################  PLOTTING THE MODELED SCENARIO ######################################


#--- Graphical functions ...........................................................................

fullHeadArrow=function(x0, y0, x1, y1, length, angle, color="black",  weight=1) {
  arrows(x0, y0, x1, y1, length, angle, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.80, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.60, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.40, code=2, lty=1, col=color, lwd=weight) 
  arrows(x0, y0, x1, y1, length, angle*0.20, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.10, code=2, lty=1, col=color, lwd=weight)
}
drawTriangle=function(growth, x, y, size, aspRatio, color) {
  if (growth>0) {
    x0=x; y0=y
    x1=x-size/2; y1=y+size/2*aspRatio
    x2=x+size/2; y2=y+size/2*aspRatio
  } else {
    x0=x-size/2; y0=y
    x1=x; y1=y+size/2*aspRatio
    x2=x+size/2; y2=y
  }
  polygon(c(x0, x1, x2, x0), c(y0, y1, y2, y0), col=color)
  return(y+size/2*aspRatio)
}


#--- Computing maximum current pop size ............................................................
maxPopSize=popSizes[1]
minPopSize=popSizes[1]
if (length(popSizes)>1){
  for (i in 2:length(popSizes)) {
    if (popSizes[i]>maxPopSize) maxPopSize=popSizes[i];
    if (popSizes[i]<minPopSize) minPopSize=popSizes[i];
  }
}

#--- Find min and max pop sizes over the whole population history ..................................
ps=popSizes
isGrowth=FALSE
if (numHistEvents) {
  #Need to keep track ofgrowth rates over time
  gRates=growthRates
  prevTime=0
  for (i in 1:numHistEvents) {
    he=histEvents[i,]
    curTime=he[1]
    sink=he[3]+1 
    resize=he[5]
    growth=he[6]
    
    #Begin by resizing pop sizes due to growth since last event
    for (j in 1:length(ps)) {
      ps[j]=ps[j]*exp(gRates[j]*(curTime-prevTime))
      if (ps[j]>maxPopSize) maxPopSize=ps[j]
      if (ps[j]<minPopSize) minPopSize=ps[j]
      if (ps[j]==0) print(paste("Deme ", j-1, " reaches size zero at time ",  curTime, " due to large negative growth (", gRates[j], ")", sep=""))
      if (is.infinite(ps[j])) print(paste("Deme ", j-1, " reaches infinite size at time ",  curTime, " due to large positive growth (", gRates[j], ")", sep=""))
    }
    prevTime=curTime
    
    #Handle transformed keep keyword
    if(growth!=-9999) {
      gRates[sink]=growth 
    }
    
    #Implement resize
    ps[sink]=ps[sink]*resize
    if (ps[sink]>maxPopSize) maxPopSize=ps[sink]
    if (ps[sink]<minPopSize) minPopSize=ps[sink]
  }
}

#-- Function to compute the circle radius for a given pop size ....................................
interpolRadius=function(curSize, minSize, maxSize, minRadius, maxRadius, logScale) {
  if(logScale) {
    minSize=log10(minSize)
    maxSize=log10(maxSize)
    curSize=log10(curSize)
  }
  curRadius=minRadius+(curSize-minSize)*(maxRadius-minRadius)/(maxSize-minSize)
  curRadius
}

#==============================                     ========================
#==============================   BEGIN MAIN PLOT   ========================
#==============================                     ========================

library("plotrix")
library("diagram")

## Loading required package: shape

## 
## Attaching package: 'diagram'

## The following object is masked from 'package:sp':
## 
##     coordinates

if (printPDF) {
  pdf(pdfFileName, width=pdf.x.size, height=pdf.y.size)
} 

par(xpd=F, mar=c(6,4,3,0.5))

maxRadius=maxRadius*(numSamples+2)
minRadius=minRadius*(numSamples+2)

title=parFileName

plot(x=1:numSamples, type="n", xlab="", ylab="time (gen)", xlim=c(-0.5, numSamples+1.5), cex.main=0.8,
     ylim=c(0, yTimeLimit), main=title, xaxt = 'n', log=logScaleAxis, cex.axis=0.9, cex.lab=0.9)
axis(side=1, labels=c("Mig Mat",0:(numSamples-1), " \nTimes"), at=0:(numSamples+1), cex.axis=0.8)
mtext("demes", side=1,line=2, cex=0.8)

w <- par("pin")[1]/diff(par("usr")[1:2])
h <- par("pin")[2]/diff(par("usr")[3:4])
aspRatio <- w/h

# aspRatio=yTimeLimit/(numSamples+2)

#--- Plotting population circles according to their current size and sampling time
for (i in 1:numSamples) {
  curRadius=interpolRadius(popSizes[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
  #   print(curRadius)
  if (i==1) topOfCircle=sampTimes[i]+curRadius*aspRatio else topOfCircle=c(topOfCircle, sampTimes[i]+curRadius*aspRatio)
  curColor=popCol
  if (sampSizes[i]==0) curColor="white"
  if (inbrCoeff[i]==0) {
    draw.circle(i, sampTimes[i], radius = curRadius, col=curColor, border=popBorderCol) 
  } else {
    floating.pie(i, sampTimes[i], c(inbrCoeff[i], 1-inbrCoeff[i]), radius = curRadius, 
                 col=c(inbreedColor,curColor), startpos=pi, border=popBorderCol)    
  }
  if (sampTimes[i]>0) {
    text(i, sampTimes[i]-curRadius*aspRatio, labels=sampTimes[i], cex=timeProp, col=oldSampColor, pos=1)
  }
  #Draw a vertical arrow in case of pop growth
  curGrowthRate=growthRates[i]
  
  #--- Draw gtrowing or shrinking triangle on top of pop circle to reflect growth type
  if (curGrowthRate!=0) {
    arLength=0.15
    topOfCircle[i]=drawTriangle(curGrowthRate, x=i, y=topOfCircle[i], arLength, aspRatio, color=growthCol)
  }
}

#--- Handle first migration matrix .......................
curMigMatNum=0
curvature=0.0075*last.he.time
text(0-migrOffset, 0, labels=0, cex=migMatNameProp, col=migrMatCol)
if (numMigMat) {
  curMigMat=migMats[1][[1]]
  for (sink in 1:numSamples) {
    for (sourc in 1:numSamples) {
      if (sink!=sourc & curMigMat[sourc, sink]>0)  {
        differ=sourc-sink
        curvedarrow(from=c(sourc, 0), to=c(sink,0), curve=curvature*(abs(differ)*0.55^abs(differ)), arr.adj=1, 
                    arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY, 
                    lcol=migrMatCol, arr.length=arrowLength)   
        if (plotMigrRates) {
          curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)        
          if (differ>0) {
            text(sink+abs(differ)/2, aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
          } else {
            text(sourc+abs(differ)/2, -aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
          } 
        }
      }
    }
  }
}

#---Draw all events on the population tree
lastTime=0
activePops=1:numSamples
numActivePops=numSamples
lastSink=-1
if (numHistEvents) {
  for (i in 1:numHistEvents) {
    #Extract historical event
    he=histEvents[i,]
    he.time=he[1]
    he.source=he[2]+1 #+1 is due to the use of base 0 for deme number in fsc
    he.sink=he[3]+1
    he.migr=he[4]
    he.resize=he[5]
    he.growth=he[6]
    if(he.growth==-9999) { #handle transformed keep keyword
      he.growth=growthRates[he.sink] #Keep current growth rate
    }
    he.migrMat=he[7]
    
    if(he.migrMat==-9999) { #handle transformed keep keyword
      he.migrMat=curMigMatNum
    }
    
    #-- Draw event time
    if (i%%2) {
      slide=timeOffset
    }
    else {
      slide=-timeOffset
    }
    
    #Draw all vertical segments ....................
    for (p in 1:length(activePops)) {
      #Handle population with older sampling times
      if (sampTimes[activePops[p]]>0) {
        minTime=max(lastTime, topOfCircle[activePops[p]]);
      } else {
        minTime=lastTime;
      }
      if (he.time>sampTimes[activePops[p]] & he.time>topOfCircle[activePops[p]]) {
        if (activePops[p]==lastSink | i==1) {
          segments(activePops[p], topOfCircle[activePops[p]], activePops[p], he.time)    
        } else {
          segments(activePops[p], minTime, activePops[p], he.time)        
        }
      }
    }     
    
    #Handle growth rate changes since last event ......................
    #Update pop sizes according to current growth rates
    if (i>1) {
      prev.he=histEvents[i-1,] 
      branchLength=he.time-lastTime
      for (p in 1:length(activePops)) {
        curPop=activePops[p]
        popSizes[curPop]=popSizes[curPop]*exp(growthRates[curPop]*branchLength);
      }
    }
    #Update growth rate
    growthRates[he.sink]=he.growth
    
    
    lastTime=he.time
    #Handle resize of sink pop ........................
    popSizes[he.sink]=popSizes[he.sink]*he.resize
    
    curRadius=interpolRadius(popSizes[he.sink], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
    #     popRadius[he.sink]=curRadius
    topOfCircle[he.sink]=he.time+curRadius*aspRatio
    draw.circle(he.sink, he.time, radius = curRadius, col=popCol, border=popBorderCol)
    
    
    #--- Draw growing or shrinking triangle on top of pop circle to reflect growth type
    if (he.growth!=0) {
      arLength=0.15
      topOfCircle[he.sink]=drawTriangle(he.growth, x=he.sink, y=topOfCircle[he.sink], arLength, aspRatio, color=growthCol)
    }
    
    #Handle population fusion .........................
    if (he.migr>=1 & he.sink!=he.source) { #This is a population fusion
      if (numActivePops==numSamples) removedPops=he.source else removedPops=c(removedPops, he.source)
      numActivePops=numActivePops-1
      activePops=(1:numSamples)[-removedPops]
      #Draw connecting arrows from source to sink
      fullHeadArrow(he.source, he.time, he.sink, he.time, length=0.15, angle=20)      
      #Redraw time with the right color
      text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=popFusionColor)
    } else {
      #--- Handle admixture event ........................
      if (he.migr>0 & he.migr<1) { 
        #Draw connecting arrows from source to sink      
        segments(he.source, he.time, he.sink, he.time, col=admixCol, lty=2)
        if (he.sink>he.source) {
          fullHeadArrow(he.sink-0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
        } else {
          fullHeadArrow(he.sink+0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
        }     
        #Redraw time with the right color
        text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=admixCol)
      }
      else text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=timeCol)
    }    
    
    
    #--- Handle migmat change  ............................
    if (i!=numHistEvents) nextTime=histEvents[i+1,1] else {
      nextTime=yTimeLimit
    }
    time2DrawArrows=(he.time+nextTime)/2
    
    if (he.migrMat!=curMigMatNum) {
      if (he.migrMat>-1) {
        curMigMat=migMats[he.migrMat+1][[1]]
        for (sink in 1:numSamples) {
          for (sourc in 1:numSamples) {
            if (sink!=sourc & curMigMat[sourc, sink]>0)  {  
              differ=sourc-sink
              curvedarrow(from=c(sourc, time2DrawArrows), to=c(sink, time2DrawArrows), curve=curvature*(abs(differ)*0.55^abs(differ)), 
                          arr.adj=1, arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY, 
                          lcol=migrMatCol, arr.length=arrowLength)  
              if (plotMigrRates) {
                #Write Nm values
                curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)
                if (differ>0) {
                  text(sink+abs(differ)/2, time2DrawArrows+aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
                } else {
                  text(sourc+abs(differ)/2, time2DrawArrows-aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
                } 
              }
            }
          }
          # curMigMatNum=he.migrMat
        }
      } 
      #--- Draw separation between migration matrices numbers on the left
      segments(-migMatLineLength/2, he.time, migMatLineLength/2, he.time, lty=3, col=migrMatCol)
    } 
    curMigMatNum=he.migrMat
    
    #Output current valid migration matrix      
    if (i%%2) {
      slide=migrOffset
    }
    else {
      slide=-migrOffset
    }   
    if (he.migrMat>-1) {
      migText=he.migrMat
      curCex=migMatNameProp
    } else {
      migText="nomig"
      curCex=migMatNameProp/2
      slide=slide*2
    }
    text(0+slide, time2DrawArrows, labels=migText, cex=curCex, col=migrMatCol)    
    
    lastSink=he.sink
  }
}

#-- Draw last branch
segments(activePops[1],topOfCircle[activePops[1]], activePops[1], yTimeLimit)

#==============================   PLOT LEGENDS IN MARGINS   ========================

#Compute space available in margin
minY.coo=grconvertY(0, from="nic", to="user")

par(xpd=NA)
#--- Draw population size scale with circles of different sizes .................
maxOrder=ceiling(log10(maxPopSize))
minOrder=floor(log10(minPopSize))
popSizeRadius=10^(maxOrder:minOrder)
winWidth=numSamples+2
ypos=3/4*minY.coo
text(x=-winWidth/10*1.2, y=ypos, labels="Pop. \nsizes ", cex=.8, pos=2)

for (i in 1:length(popSizeRadius)) {
  curRadius=interpolRadius(popSizeRadius[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
  #   print(curRadius)
  if (curRadius>0) {
    xpos=-winWidth/10+(i-1)*winWidth/10
    draw.circle(xpos, ypos, radius=curRadius, col=popCol, border=popBorderCol)
  }
  text(xpos, ypos-abs(ypos)*0.1, popSizeRadius[i], cex=0.7, pos=1, col="black")
}

#--- Legend for growing or shrinking populations ...............................
if (isGrowth) {
  x=winWidth-1.5*winWidth/10; y=ypos+abs(ypos)*0.2
  text(x, y-abs(ypos)*0.1, labels="Populations", cex=0.8)
  x=winWidth-2*winWidth/10; y=ypos-abs(ypos)*0.1
  drawTriangle(1, x, y, size=0.15, aspRatio, color=growthCol)
  text(x, y-abs(ypos)*0.1, labels="growing", pos=NULL, cex=0.7)
  x=winWidth-winWidth/10
  drawTriangle(-1, x, y, size=0.15, aspRatio, color=growthCol)
  text(x, y-abs(ypos)*0.1, labels="shrinking", pos=NULL, cex=0.7)  
}


#if (printPDF) dev.off()
dev.off()

## null device 
##           1

Varroa jacobsoni demographic history

The same process described for V. destructor demographic history inferences was independently applied for its sister species V. jacobsoni.

Likelihood distribution for the four VJ SFS datasets

Contrary to V. destructor for which the scenario 6 was undeniably the most likely according to both AIC and model likelihood distribution, both scenarios 3 and 6 seems likely in the case of recent host switch by V. jacobsoni.

As a consequence, we repeated the same process of bootstraps computations as done for V. destructor but for both scenarios: DIV + BOT and IM + BOT.

vj.lhoods1050 <- read.csv("fsc_run/lhoods_1050SNPS_vjpng.csv", header = TRUE)
vj.lhoods10500 <- read.csv("fsc_run/lhoods_10500SNPS_vjpng.csv", header = TRUE)
vj.lhoods35000 <- read.csv("fsc_run/lhoods_35000SNPS_vjpng.csv", header = TRUE)
vj.lhoods105000 <- read.csv("fsc_run/lhoods_105000SNPS_vjpng.csv", header = TRUE)

par(mfrow=c(1,4))
#Plot 1 for likelihoods with smaller subset
boxplot(range = 0,vj.lhoods1050$scenario1_div [3:102], vj.lhoods1050$scenario2_divGwt2n[3:102], vj.lhoods1050$scenario3_divGwt[3:102], vj.lhoods1050$scenario4_IM[3:102], vj.lhoods1050$scenario5_IMGwt2n[3:102], vj.lhoods1050$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vj.lhoods1050$scenario1_div [2],col = "blue", pch = 20)
points(2,vj.lhoods1050$scenario2_divGwt2n [2],col = "blue", pch = 20)
points(3,vj.lhoods1050$scenario3_divGwt [2],col = "blue", pch = 20)
points(4,vj.lhoods1050$scenario4_IM [2],col = "blue", pch = 20)
points(5,vj.lhoods1050$scenario5_IMGwt2n [2],col = "blue", pch = 20)
points(6,vj.lhoods1050$scenario6_IMGwt [2],col = "blue", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. jacobsoni 472 SNPS")

#Plot 2 for likelihoods with second smaller subset
boxplot(range = 0,vj.lhoods10500$scenario1_div [3:102], vj.lhoods10500$scenario2_divGwt2n[3:102], vj.lhoods10500$scenario3_divGwt[3:102], vj.lhoods10500$scenario4_IM[3:102], vj.lhoods10500$scenario5_IMGwt2n[3:102], vj.lhoods10500$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vj.lhoods10500$scenario1_div [2],col = "blue", pch = 20)
points(2,vj.lhoods10500$scenario2_divGwt2n [2],col = "blue", pch = 20)
points(3,vj.lhoods10500$scenario3_divGwt [2],col = "blue", pch = 20)
points(4,vj.lhoods10500$scenario4_IM [2],col = "blue", pch = 20)
points(5,vj.lhoods10500$scenario5_IMGwt2n [2],col = "blue", pch = 20)
points(6,vj.lhoods10500$scenario6_IMGwt [2],col = "blue", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. jacobsoni 4,713 SNPS")

#Plot 3 for likelihoods with third subset
boxplot(range = 0,vj.lhoods35000$scenario1_div [3:102], vj.lhoods35000$scenario2_divGwt2n[3:102], vj.lhoods35000$scenario3_divGwt[3:102], vj.lhoods35000$scenario4_IM[3:102], vj.lhoods35000$scenario5_IMGwt2n[3:102], vj.lhoods35000$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vj.lhoods35000$scenario1_div [2],col = "blue", pch = 20)
points(2,vj.lhoods35000$scenario2_divGwt2n [2],col = "blue", pch = 20)
points(3,vj.lhoods35000$scenario3_divGwt [2],col = "blue", pch = 20)
points(4,vj.lhoods35000$scenario4_IM [2],col = "blue", pch = 20)
points(5,vj.lhoods35000$scenario5_IMGwt2n [2],col = "blue", pch = 20)
points(6,vj.lhoods35000$scenario6_IMGwt [2],col = "blue", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. jacobsoni 15,363 SNPS")

#Plot 4 for likelihoods with larger subset
boxplot(range = 0,vj.lhoods105000$scenario1_div [3:102], vj.lhoods105000$scenario2_divGwt2n[3:102], vj.lhoods105000$scenario3_divGwt[3:102], vj.lhoods105000$scenario4_IM[3:102], vj.lhoods105000$scenario5_IMGwt2n[3:102], vj.lhoods105000$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vj.lhoods105000$scenario1_div [2],col = "blue", pch = 20)
points(2,vj.lhoods105000$scenario2_divGwt2n [2],col = "blue", pch = 20)
points(3,vj.lhoods105000$scenario3_divGwt [2],col = "blue", pch = 20)
points(4,vj.lhoods105000$scenario4_IM [2],col = "blue", pch = 20)
points(5,vj.lhoods105000$scenario5_IMGwt2n [2],col = "blue", pch = 20)
points(6,vj.lhoods105000$scenario6_IMGwt [2],col = "blue", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. jacobsoni 46,416 SNPS")

Scenario 3: DIV + BOT bootstraps computation

# import the file containing the 100 best runs
vj.boot3 <- read.csv("fsc_run/png_105000nopv_scen3_BOOT.csv", header = TRUE)

# Effective size for Varroa mites on original host Apis cerana
N_VjAc1.boot = boot(vj.boot3$NVJAC1, my.mean, 10000)
my.mean(vj.boot3$NVJAC1,1:length(vj.boot3$NVJAC1))

## [1] 228.34

boot.ci(N_VjAc1.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = N_VjAc1.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (208.7, 248.0 )   (208.0, 247.6 )  
## 
## Level     Percentile            BCa          
## 95%   (209.1, 248.6 )   (209.5, 249.0 )  
## Calculations and Intervals on Original Scale

# Effective size for Varroa mites on new host Apis mellifera
N_VjAm0.boot = boot(vj.boot3$NVJAM0, my.mean, 10000)
my.mean(vj.boot3$NVJAM0,1:length(vj.boot3$NVJAM0))

## [1] 6708.25

boot.ci(N_VjAm0.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = N_VjAm0.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (6231, 7193 )   (6229, 7191 )  
## 
## Level     Percentile            BCa          
## 95%   (6226, 7188 )   (6235, 7198 )  
## Calculations and Intervals on Original Scale

# Number of generations when new population of Varroa mites was founded
TBOT.boot = boot(vj.boot3$TBOT, my.mean, 10000)
my.mean(vj.boot3$TBOT,1:length(vj.boot3$TBOT))

## [1] 58.66

boot.ci(TBOT.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = TBOT.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (53.44, 63.74 )   (53.34, 63.67 )  
## 
## Level     Percentile            BCa          
## 95%   (53.65, 63.98 )   (53.67, 64.01 )  
## Calculations and Intervals on Original Scale

# Growth rate of Varroa mites on new host Apis mellifera
GAM.boot = boot(vj.boot3$GAM, my.mean, 10000)
my.mean(vj.boot3$GAM,1:length(vj.boot3$GAM))

## [1] -0.156716

boot.ci(GAM.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = GAM.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (-0.1730, -0.1405 )   (-0.1722, -0.1400 )  
## 
## Level     Percentile            BCa          
## 95%   (-0.1734, -0.1413 )   (-0.1750, -0.1423 )  
## Calculations and Intervals on Original Scale

# Effective size for Varroa mites able to jump on A. mellifera to create new population
JUMP.boot = boot(vj.boot3$JUMP, my.mean, 10000)
my.mean(vj.boot3$JUMP,1:length(vj.boot3$JUMP))

## [1] 348.86

boot.ci(JUMP.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = JUMP.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (313.9, 384.3 )   (313.2, 382.8 )  
## 
## Level     Percentile            BCa          
## 95%   (314.9, 384.5 )   (317.2, 388.7 )  
## Calculations and Intervals on Original Scale

Draw parameters for scenario 3 with mean bootstraps

#...........................................................................................
# (c) Laurent Excoffier and Vitor Sousa June-November 2015
# 
# Small R program to draw the evolutonary scenario described by a given par file
# This is mainly for visual checking that the modeled scenario corresponds to   
# what was intended
#
#...........................................................................................
#
args=commandArgs(TRUE)
print(args)

## character(0)

#Expects par file name on command line

if(length(args)) {
  parFileName=args[1]
} else {  
  #REPLACE BY THE NAME OF THE PAR FILE YOU WANT TO ANALYSE
  parFileName="fsc_run/VJM3_downsample_divGwt_29_boot.par"
}

separator=" "
migrMatCol="coral"
admixCol="blue"
popFusionColor="black"
popCol="lightgrey"
popBorderCol="black"
inbreedColor="lightblue"
oldSampColor="olivedrab4"
timeCol="tan4"
growthCol="hotpink3"
propLastsegment=0.05
migMatNameProp=0.8
migMatLineLength=0.3
timeProp=0.6
maxRadius=1/40
minRadius=maxRadius/3
arrowLength=0.2
logScaleAxis=""
timeOffset=0.25
migrOffset=0.1
curvedArrowLTY=1
drawLogPopSize=T
plotMigrRates=T
migrRateTextSizeCEX=0.5

#Define plot area size for PDF
pdf.x.size=7
pdf.y.size=10

rescalingFactor=1.0 #Don't touch this!


printPDF=F

# parFile=readLines(con=parFileName)

###########################      Reading par file      #########################

parFile=scan(parFileName, character(0), sep = "\n", strip.white = TRUE) # separate each line

pdfFileName=paste(parFileName, ".pdf", sep="")

#--- Clean input file by removing consecutive separators, and keep ---------------

#--- Function to remove separators within a string
removeTrailingSep=function(string, sep) {
  temp=strsplit(string, split=sep)
  temp2=temp[[1]][nchar(temp[[1]])>0]
  cleanStr=temp2[1]
  if (length(temp2)>1) {
    for (i in 2:length(temp2)) {
      cleanStr=paste(cleanStr, temp2[i], sep=sep)
    }
  }
  cleanStr
}

#--- Replace Keep by -9999
replaceKeep=function(string) {
  if (grepl("keep", string)) {
    return(gsub("keep", "-9999", string))
  }
  return(string)
}

#Remove both multiple consecutive whitespace and tabs and replace keep keyword
for (i in 1:length(parFile)) {
  parFile[i]=removeTrailingSep(parFile[i], sep='\t')
  parFile[i]=removeTrailingSep(parFile[i], sep=' ')
  parFile[i]=replaceKeep(parFile[i])
}

#-------------------------------------------------------------------------------

#--- Get number of samples on line 2 -----
l.numsamples=parFile[2]
sp.l=unlist(strsplit(l.numsamples, split=' '))
tab.l=unlist(strsplit(l.numsamples, split='\t'))
if (length(sp.l)>=length(tab.l)) {
  numSamples=as.numeric(sp.l[1])
  separator=" "
} else {
  numSamples=as.numeric(tab.l[1])
  separator="\t"
}

#--- Reading numbers on separate lines -----
getNumbers=function(start, parFile, numSamples) {  
  for (i in 1:numSamples) {
    curnum=as.numeric(unlist(strsplit(parFile[start+i], split=separator))[1])
    if (i==1) {
      num=curnum 
    } else {
      num=c(num, curnum)
    }
  }
  num
}

#--- Get population sizes -----------------
start=3
popSizes=getNumbers(start, parFile, numSamples)
#Rescaling pop sizes
popSizes=round(popSizes*rescalingFactor, digits=0)

iniPopSizes=popSizes

#--- Get sample sizes -----------------

readSampleSizesTimesAndInbreedingLevel=function(start, parFile, numSamples) {
  for (i in 1:numSamples) {
    curLine=unlist(strsplit(parFile[start+i], split=separator))
    curSampSize=as.numeric(curLine[1])
    curSampTime=0
    curInbreeding=0
    if (length(curLine)>1) curSampTime=as.numeric(curLine[2])
    if (length(curLine)>2) curInbreeding=as.numeric(curLine[3])
    if (i==1) {
      sampSize=curSampSize
      sampTime=curSampTime
      inbreeding=curInbreeding
    } else {
      sampSize=c(sampSize,curSampSize)
      sampTime=c(sampTime,curSampTime)
      inbreeding=c(inbreeding,curInbreeding)
    }
  }
  list(ss=sampSize, st=sampTime, inb=inbreeding)
}

start=start+numSamples+1
# sampSizes=getNumbers(start, parFile, numSamples)
sampSizesStats=readSampleSizesTimesAndInbreedingLevel(start, parFile, numSamples)
#Rescaling sample times
sampSizesStats$st=round(sampSizesStats$st*rescalingFactor, digits=0)

sampSizes=sampSizesStats$ss
sampTimes=sampSizesStats$st
inbrCoeff=sampSizesStats$inb

#--- Get growth rates -----------------
start=start+numSamples+1
growthRatesInitial=getNumbers(start, parFile, numSamples)
# save this into growthRates which will be used and updated when printing historical events
growthRates=growthRatesInitial

#--- Get number of migration matrices -----------------
start=start+numSamples+1
numMigMat=as.numeric(unlist(strsplit(parFile[start+1], split=separator))[1])

#--- Read migration matrix
readMigMat=function(start, parFile, numSamples) {
  for (i in 1:numSamples) {
    curmigs=as.numeric(unlist(strsplit(parFile[start+i], split=separator)))
    if (i==1) {
      migs=curmigs 
    } else {
      migs=rbind(migs, curmigs)
    }
  }
  rownames(migs)=1:numSamples
  migs 
}

#--- Get migration matrices as a list --------------
start=start+2
migMats=list()
if (numMigMat>0) {
  for (i in 1:numMigMat) {  
    curMigMat=readMigMat(start, parFile, numSamples) 
    migMats[[i]]=curMigMat
    start=start+numSamples+1
  }
}

#Rescaling migration rates
if (numMigMat>0) {
  for (i in 1:numMigMat) {
    migMats[[i]]=migMats[[i]]/rescalingFactor
  }
}

#--- Get number of historical events
start=start+1
numHistEvents=as.numeric(unlist(strsplit(parFile[start], split=separator))[1])

###################### HISTORICAL EVENTS HANDLING ##############################

#..... Read Historical Event .......
last.he.time=0
if (numHistEvents>0) {
  for (i in 1:numHistEvents) {
    start=start+1
    #Take care of nomig keyword
    nomig=F
    if (grepl("nomig", parFile[start])) {
      nomig=T
      gsub("nomig", "", parFile[start])
    }
    curHE=as.numeric(unlist(strsplit(parFile[start], split=separator)))
    if (nomig) curHE[7]=-1 
    #Rescaling time of event
    curHE[1]=round(curHE[1]*rescalingFactor, digits=0)
    if (i==1) {
      histEvents=curHE
      last.he.time=curHE[1]
    } else {
      histEvents=rbind(histEvents, curHE)
      if (histEvents[i,1] > last.he.time) last.he.time=histEvents[i,1] 
    }
  }
  if (numHistEvents>1) rownames(histEvents)=1:numHistEvents
}
names(last.he.time)=""
last.he.time=as.numeric(last.he.time)

yTimeLimit=0
if (last.he.time!=0) yTimeLimit=last.he.time*(1+propLastsegment) #Add propLastsegment to y axis after last event (to draw stuff)

#Reorder events by their times
if (numHistEvents>1) histEvents=histEvents[order(histEvents[,1],decreasing=FALSE),] else {
  if (numHistEvents==1) histEvents=matrix(histEvents, nrow=1,  byrow = T)
}

endReadParFile=start



##############################  PLOTTING THE MODELED SCENARIO ######################################


#--- Graphical functions ...........................................................................

fullHeadArrow=function(x0, y0, x1, y1, length, angle, color="black",  weight=1) {
  arrows(x0, y0, x1, y1, length, angle, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.80, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.60, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.40, code=2, lty=1, col=color, lwd=weight) 
  arrows(x0, y0, x1, y1, length, angle*0.20, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.10, code=2, lty=1, col=color, lwd=weight)
}
drawTriangle=function(growth, x, y, size, aspRatio, color) {
  if (growth>0) {
    x0=x; y0=y
    x1=x-size/2; y1=y+size/2*aspRatio
    x2=x+size/2; y2=y+size/2*aspRatio
  } else {
    x0=x-size/2; y0=y
    x1=x; y1=y+size/2*aspRatio
    x2=x+size/2; y2=y
  }
  polygon(c(x0, x1, x2, x0), c(y0, y1, y2, y0), col=color)
  return(y+size/2*aspRatio)
}


#--- Computing maximum current pop size ............................................................
maxPopSize=popSizes[1]
minPopSize=popSizes[1]
if (length(popSizes)>1){
  for (i in 2:length(popSizes)) {
    if (popSizes[i]>maxPopSize) maxPopSize=popSizes[i];
    if (popSizes[i]<minPopSize) minPopSize=popSizes[i];
  }
}

#--- Find min and max pop sizes over the whole population history ..................................
ps=popSizes
isGrowth=FALSE
if (numHistEvents) {
  #Need to keep track ofgrowth rates over time
  gRates=growthRates
  prevTime=0
  for (i in 1:numHistEvents) {
    he=histEvents[i,]
    curTime=he[1]
    sink=he[3]+1 
    resize=he[5]
    growth=he[6]
    
    #Begin by resizing pop sizes due to growth since last event
    for (j in 1:length(ps)) {
      ps[j]=ps[j]*exp(gRates[j]*(curTime-prevTime))
      if (ps[j]>maxPopSize) maxPopSize=ps[j]
      if (ps[j]<minPopSize) minPopSize=ps[j]
      if (ps[j]==0) print(paste("Deme ", j-1, " reaches size zero at time ",  curTime, " due to large negative growth (", gRates[j], ")", sep=""))
      if (is.infinite(ps[j])) print(paste("Deme ", j-1, " reaches infinite size at time ",  curTime, " due to large positive growth (", gRates[j], ")", sep=""))
    }
    prevTime=curTime
    
    #Handle transformed keep keyword
    if(growth!=-9999) {
      gRates[sink]=growth 
    }
    
    #Implement resize
    ps[sink]=ps[sink]*resize
    if (ps[sink]>maxPopSize) maxPopSize=ps[sink]
    if (ps[sink]<minPopSize) minPopSize=ps[sink]
  }
}

#-- Function to compute the circle radius for a given pop size ....................................
interpolRadius=function(curSize, minSize, maxSize, minRadius, maxRadius, logScale) {
  if(logScale) {
    minSize=log10(minSize)
    maxSize=log10(maxSize)
    curSize=log10(curSize)
  }
  curRadius=minRadius+(curSize-minSize)*(maxRadius-minRadius)/(maxSize-minSize)
  curRadius
}

#==============================                     ========================
#==============================   BEGIN MAIN PLOT   ========================
#==============================                     ========================

library("plotrix")
library("diagram")

if (printPDF) {
  pdf(pdfFileName, width=pdf.x.size, height=pdf.y.size)
} 

par(xpd=F, mar=c(6,4,3,0.5))

maxRadius=maxRadius*(numSamples+2)
minRadius=minRadius*(numSamples+2)

title=parFileName

plot(x=1:numSamples, type="n", xlab="", ylab="time (gen)", xlim=c(-0.5, numSamples+1.5), cex.main=0.8,
     ylim=c(0, yTimeLimit), main=title, xaxt = 'n', log=logScaleAxis, cex.axis=0.9, cex.lab=0.9)
axis(side=1, labels=c("Mig Mat",0:(numSamples-1), " \nTimes"), at=0:(numSamples+1), cex.axis=0.8)
mtext("demes", side=1,line=2, cex=0.8)

w <- par("pin")[1]/diff(par("usr")[1:2])
h <- par("pin")[2]/diff(par("usr")[3:4])
aspRatio <- w/h

# aspRatio=yTimeLimit/(numSamples+2)

#--- Plotting population circles according to their current size and sampling time
for (i in 1:numSamples) {
  curRadius=interpolRadius(popSizes[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
  #   print(curRadius)
  if (i==1) topOfCircle=sampTimes[i]+curRadius*aspRatio else topOfCircle=c(topOfCircle, sampTimes[i]+curRadius*aspRatio)
  curColor=popCol
  if (sampSizes[i]==0) curColor="white"
  if (inbrCoeff[i]==0) {
    draw.circle(i, sampTimes[i], radius = curRadius, col=curColor, border=popBorderCol) 
  } else {
    floating.pie(i, sampTimes[i], c(inbrCoeff[i], 1-inbrCoeff[i]), radius = curRadius, 
                 col=c(inbreedColor,curColor), startpos=pi, border=popBorderCol)    
  }
  if (sampTimes[i]>0) {
    text(i, sampTimes[i]-curRadius*aspRatio, labels=sampTimes[i], cex=timeProp, col=oldSampColor, pos=1)
  }
  #Draw a vertical arrow in case of pop growth
  curGrowthRate=growthRates[i]
  
  #--- Draw gtrowing or shrinking triangle on top of pop circle to reflect growth type
  if (curGrowthRate!=0) {
    arLength=0.15
    topOfCircle[i]=drawTriangle(curGrowthRate, x=i, y=topOfCircle[i], arLength, aspRatio, color=growthCol)
  }
}

#--- Handle first migration matrix .......................
curMigMatNum=0
curvature=0.0075*last.he.time
text(0-migrOffset, 0, labels=0, cex=migMatNameProp, col=migrMatCol)
if (numMigMat) {
  curMigMat=migMats[1][[1]]
  for (sink in 1:numSamples) {
    for (sourc in 1:numSamples) {
      if (sink!=sourc & curMigMat[sourc, sink]>0)  {
        differ=sourc-sink
        curvedarrow(from=c(sourc, 0), to=c(sink,0), curve=curvature*(abs(differ)*0.55^abs(differ)), arr.adj=1, 
                    arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY, 
                    lcol=migrMatCol, arr.length=arrowLength)   
        if (plotMigrRates) {
          curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)        
          if (differ>0) {
            text(sink+abs(differ)/2, aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
          } else {
            text(sourc+abs(differ)/2, -aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
          } 
        }
      }
    }
  }
}

#---Draw all events on the population tree
lastTime=0
activePops=1:numSamples
numActivePops=numSamples
lastSink=-1
if (numHistEvents) {
  for (i in 1:numHistEvents) {
    #Extract historical event
    he=histEvents[i,]
    he.time=he[1]
    he.source=he[2]+1 #+1 is due to the use of base 0 for deme number in fsc
    he.sink=he[3]+1
    he.migr=he[4]
    he.resize=he[5]
    he.growth=he[6]
    if(he.growth==-9999) { #handle transformed keep keyword
      he.growth=growthRates[he.sink] #Keep current growth rate
    }
    he.migrMat=he[7]
    
    if(he.migrMat==-9999) { #handle transformed keep keyword
      he.migrMat=curMigMatNum
    }
    
    #-- Draw event time
    if (i%%2) {
      slide=timeOffset
    }
    else {
      slide=-timeOffset
    }
    
    #Draw all vertical segments ....................
    for (p in 1:length(activePops)) {
      #Handle population with older sampling times
      if (sampTimes[activePops[p]]>0) {
        minTime=max(lastTime, topOfCircle[activePops[p]]);
      } else {
        minTime=lastTime;
      }
      if (he.time>sampTimes[activePops[p]] & he.time>topOfCircle[activePops[p]]) {
        if (activePops[p]==lastSink | i==1) {
          segments(activePops[p], topOfCircle[activePops[p]], activePops[p], he.time)    
        } else {
          segments(activePops[p], minTime, activePops[p], he.time)        
        }
      }
    }     
    
    #Handle growth rate changes since last event ......................
    #Update pop sizes according to current growth rates
    if (i>1) {
      prev.he=histEvents[i-1,] 
      branchLength=he.time-lastTime
      for (p in 1:length(activePops)) {
        curPop=activePops[p]
        popSizes[curPop]=popSizes[curPop]*exp(growthRates[curPop]*branchLength);
      }
    }
    #Update growth rate
    growthRates[he.sink]=he.growth
    
    
    lastTime=he.time
    #Handle resize of sink pop ........................
    popSizes[he.sink]=popSizes[he.sink]*he.resize
    
    curRadius=interpolRadius(popSizes[he.sink], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
    #     popRadius[he.sink]=curRadius
    topOfCircle[he.sink]=he.time+curRadius*aspRatio
    draw.circle(he.sink, he.time, radius = curRadius, col=popCol, border=popBorderCol)
    
    
    #--- Draw growing or shrinking triangle on top of pop circle to reflect growth type
    if (he.growth!=0) {
      arLength=0.15
      topOfCircle[he.sink]=drawTriangle(he.growth, x=he.sink, y=topOfCircle[he.sink], arLength, aspRatio, color=growthCol)
    }
    
    #Handle population fusion .........................
    if (he.migr>=1 & he.sink!=he.source) { #This is a population fusion
      if (numActivePops==numSamples) removedPops=he.source else removedPops=c(removedPops, he.source)
      numActivePops=numActivePops-1
      activePops=(1:numSamples)[-removedPops]
      #Draw connecting arrows from source to sink
      fullHeadArrow(he.source, he.time, he.sink, he.time, length=0.15, angle=20)      
      #Redraw time with the right color
      text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=popFusionColor)
    } else {
      #--- Handle admixture event ........................
      if (he.migr>0 & he.migr<1) { 
        #Draw connecting arrows from source to sink      
        segments(he.source, he.time, he.sink, he.time, col=admixCol, lty=2)
        if (he.sink>he.source) {
          fullHeadArrow(he.sink-0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
        } else {
          fullHeadArrow(he.sink+0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
        }     
        #Redraw time with the right color
        text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=admixCol)
      }
      else text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=timeCol)
    }    
    
    
    #--- Handle migmat change  ............................
    if (i!=numHistEvents) nextTime=histEvents[i+1,1] else {
      nextTime=yTimeLimit
    }
    time2DrawArrows=(he.time+nextTime)/2
    
    if (he.migrMat!=curMigMatNum) {
      if (he.migrMat>-1) {
        curMigMat=migMats[he.migrMat+1][[1]]
        for (sink in 1:numSamples) {
          for (sourc in 1:numSamples) {
            if (sink!=sourc & curMigMat[sourc, sink]>0)  {  
              differ=sourc-sink
              curvedarrow(from=c(sourc, time2DrawArrows), to=c(sink, time2DrawArrows), curve=curvature*(abs(differ)*0.55^abs(differ)), 
                          arr.adj=1, arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY, 
                          lcol=migrMatCol, arr.length=arrowLength)  
              if (plotMigrRates) {
                #Write Nm values
                curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)
                if (differ>0) {
                  text(sink+abs(differ)/2, time2DrawArrows+aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
                } else {
                  text(sourc+abs(differ)/2, time2DrawArrows-aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
                } 
              }
            }
          }
          # curMigMatNum=he.migrMat
        }
      } 
      #--- Draw separation between migration matrices numbers on the left
      segments(-migMatLineLength/2, he.time, migMatLineLength/2, he.time, lty=3, col=migrMatCol)
    } 
    curMigMatNum=he.migrMat
    
    #Output current valid migration matrix      
    if (i%%2) {
      slide=migrOffset
    }
    else {
      slide=-migrOffset
    }   
    if (he.migrMat>-1) {
      migText=he.migrMat
      curCex=migMatNameProp
    } else {
      migText="nomig"
      curCex=migMatNameProp/2
      slide=slide*2
    }
    text(0+slide, time2DrawArrows, labels=migText, cex=curCex, col=migrMatCol)    
    
    lastSink=he.sink
  }
}

#-- Draw last branch
segments(activePops[1],topOfCircle[activePops[1]], activePops[1], yTimeLimit)

#==============================   PLOT LEGENDS IN MARGINS   ========================

#Compute space available in margin
minY.coo=grconvertY(0, from="nic", to="user")

par(xpd=NA)
#--- Draw population size scale with circles of different sizes .................
maxOrder=ceiling(log10(maxPopSize))
minOrder=floor(log10(minPopSize))
popSizeRadius=10^(maxOrder:minOrder)
winWidth=numSamples+2
ypos=3/4*minY.coo
text(x=-winWidth/10*1.2, y=ypos, labels="Pop. \nsizes ", cex=.8, pos=2)

for (i in 1:length(popSizeRadius)) {
  curRadius=interpolRadius(popSizeRadius[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
  #   print(curRadius)
  if (curRadius>0) {
    xpos=-winWidth/10+(i-1)*winWidth/10
    draw.circle(xpos, ypos, radius=curRadius, col=popCol, border=popBorderCol)
  }
  text(xpos, ypos-abs(ypos)*0.1, popSizeRadius[i], cex=0.7, pos=1, col="black")
}

#--- Legend for growing or shrinking populations ...............................
if (isGrowth) {
  x=winWidth-1.5*winWidth/10; y=ypos+abs(ypos)*0.2
  text(x, y-abs(ypos)*0.1, labels="Populations", cex=0.8)
  x=winWidth-2*winWidth/10; y=ypos-abs(ypos)*0.1
  drawTriangle(1, x, y, size=0.15, aspRatio, color=growthCol)
  text(x, y-abs(ypos)*0.1, labels="growing", pos=NULL, cex=0.7)
  x=winWidth-winWidth/10
  drawTriangle(-1, x, y, size=0.15, aspRatio, color=growthCol)
  text(x, y-abs(ypos)*0.1, labels="shrinking", pos=NULL, cex=0.7)  
}


#if (printPDF) dev.off()
dev.off()

## null device 
##           1

Scenario 6: IM + BOT bootstraps computation

# import the file containing the 100 best runs
vj.boot6 <- read.csv("fsc_run/png_105000nopv_scen6_BOOT.csv", header = TRUE)

# Effective size for Varroa mites on original host Apis cerana
N_VjAc1.boot = boot(vj.boot6$NVJAC1, my.mean, 10000)
my.mean(vj.boot6$NVJAC1,1:length(vj.boot6$NVJAC1))

## [1] 150.5

boot.ci(N_VjAc1.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = N_VjAc1.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (137.2, 163.8 )   (137.0, 163.4 )  
## 
## Level     Percentile            BCa          
## 95%   (137.6, 164.0 )   (138.5, 165.0 )  
## Calculations and Intervals on Original Scale

# Effective size for Varroa mites on new host Apis mellifera
N_VjAm0.boot = boot(vj.boot6$NVJAM0, my.mean, 10000)
my.mean(vj.boot6$NVJAM0,1:length(vj.boot6$NVJAM0))

## [1] 7922.38

boot.ci(N_VjAm0.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = N_VjAm0.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (7406, 8435 )   (7398, 8431 )  
## 
## Level     Percentile            BCa          
## 95%   (7414, 8447 )   (7427, 8459 )  
## Calculations and Intervals on Original Scale

# Number of generations when new population of Varroa mites was founded
TBOT.boot = boot(vj.boot6$TBOT, my.mean, 10000)
my.mean(vj.boot6$TBOT,1:length(vj.boot6$TBOT))

## [1] 38.99

boot.ci(TBOT.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = TBOT.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (35.52, 42.43 )   (35.41, 42.28 )  
## 
## Level     Percentile            BCa          
## 95%   (35.70, 42.57 )   (35.94, 42.85 )  
## Calculations and Intervals on Original Scale

# Growth rate of Varroa mites on new host Apis mellifera
GAM.boot = boot(vj.boot6$GAM, my.mean, 10000)
my.mean(vj.boot6$GAM,1:length(vj.boot6$GAM))

## [1] -0.253887

boot.ci(GAM.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = GAM.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (-0.2776, -0.2302 )   (-0.2770, -0.2292 )  
## 
## Level     Percentile            BCa          
## 95%   (-0.2786, -0.2308 )   (-0.2799, -0.2320 )  
## Calculations and Intervals on Original Scale

# Effective size for Varroa mites able to jump on A. mellifera to create new population
JUMP.boot = boot(vj.boot6$JUMP, my.mean, 10000)
my.mean(vj.boot6$JUMP,1:length(vj.boot6$JUMP))

## [1] 295.06

boot.ci(JUMP.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = JUMP.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (268.8, 321.3 )   (268.8, 320.3 )  
## 
## Level     Percentile            BCa          
## 95%   (269.8, 321.4 )   (271.3, 323.0 )  
## Calculations and Intervals on Original Scale

N0M1.boot = boot(vj.boot6$N0M1, my.mean, 10000)
my.mean(vj.boot6$N0M1,1:length(vj.boot6$N0M1))

## [1] 0.02610133

boot.ci(N0M1.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = N0M1.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0.0166,  0.0354 )   ( 0.0160,  0.0348 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0.0174,  0.0362 )   ( 0.0183,  0.0376 )  
## Calculations and Intervals on Original Scale

N1M0.boot = boot(vj.boot6$N1M0, my.mean, 10000)
my.mean(vj.boot6$N1M0,1:length(vj.boot6$N1M0))

## [1] 0.00691537

boot.ci(N1M0.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = N1M0.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0.0044,  0.0094 )   ( 0.0044,  0.0093 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0.0045,  0.0095 )   ( 0.0048,  0.0098 )  
## Calculations and Intervals on Original Scale

MACtoAM.boot = boot(vj.boot6$MACtoAM, my.mean, 10000)
my.mean(vj.boot6$MACtoAM,1:length(vj.boot6$MACtoAM))

## [1] 5.264404e-05

boot.ci(MACtoAM.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = MACtoAM.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0.0000,  0.0001 )   ( 0.0000,  0.0001 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0.0000,  0.0001 )   ( 0.0000,  0.0001 )  
## Calculations and Intervals on Original Scale

MAMtoAC.boot = boot(vj.boot6$MAMtoAC, my.mean, 10000)
my.mean(vj.boot6$MAMtoAC,1:length(vj.boot6$MAMtoAC))

## [1] 3.453322e-06

boot.ci(MAMtoAC.boot)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = MAMtoAC.boot)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0,  0 )   ( 0,  0 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0,  0 )   ( 0,  0 )  
## Calculations and Intervals on Original Scale

Draw parameters for scenario 6 with mean bootstraps

#...........................................................................................
# (c) Laurent Excoffier and Vitor Sousa June-November 2015
# 
# Small R program to draw the evolutonary scenario described by a given par file
# This is mainly for visual checking that the modeled scenario corresponds to   
# what was intended
#
#...........................................................................................
#
args=commandArgs(TRUE)
print(args)

## character(0)

#Expects par file name on command line

if(length(args)) {
  parFileName=args[1]
} else {  
  #REPLACE BY THE NAME OF THE PAR FILE YOU WANT TO ANALYSE
  parFileName="fsc_run/VJM6_downsample_IMGwt_32_boot.par"
}

separator=" "
migrMatCol="coral"
admixCol="blue"
popFusionColor="black"
popCol="lightgrey"
popBorderCol="black"
inbreedColor="lightblue"
oldSampColor="olivedrab4"
timeCol="tan4"
growthCol="hotpink3"
propLastsegment=0.05
migMatNameProp=0.8
migMatLineLength=0.3
timeProp=0.6
maxRadius=1/40
minRadius=maxRadius/3
arrowLength=0.2
logScaleAxis=""
timeOffset=0.25
migrOffset=0.1
curvedArrowLTY=1
drawLogPopSize=T
plotMigrRates=T
migrRateTextSizeCEX=0.5

#Define plot area size for PDF
pdf.x.size=7
pdf.y.size=10

rescalingFactor=1.0 #Don't touch this!


printPDF=F

# parFile=readLines(con=parFileName)

###########################      Reading par file      #########################

parFile=scan(parFileName, character(0), sep = "\n", strip.white = TRUE) # separate each line

pdfFileName=paste(parFileName, ".pdf", sep="")

#--- Clean input file by removing consecutive separators, and keep ---------------

#--- Function to remove separators within a string
removeTrailingSep=function(string, sep) {
  temp=strsplit(string, split=sep)
  temp2=temp[[1]][nchar(temp[[1]])>0]
  cleanStr=temp2[1]
  if (length(temp2)>1) {
    for (i in 2:length(temp2)) {
      cleanStr=paste(cleanStr, temp2[i], sep=sep)
    }
  }
  cleanStr
}

#--- Replace Keep by -9999
replaceKeep=function(string) {
  if (grepl("keep", string)) {
    return(gsub("keep", "-9999", string))
  }
  return(string)
}

#Remove both multiple consecutive whitespace and tabs and replace keep keyword
for (i in 1:length(parFile)) {
  parFile[i]=removeTrailingSep(parFile[i], sep='\t')
  parFile[i]=removeTrailingSep(parFile[i], sep=' ')
  parFile[i]=replaceKeep(parFile[i])
}

#-------------------------------------------------------------------------------

#--- Get number of samples on line 2 -----
l.numsamples=parFile[2]
sp.l=unlist(strsplit(l.numsamples, split=' '))
tab.l=unlist(strsplit(l.numsamples, split='\t'))
if (length(sp.l)>=length(tab.l)) {
  numSamples=as.numeric(sp.l[1])
  separator=" "
} else {
  numSamples=as.numeric(tab.l[1])
  separator="\t"
}

#--- Reading numbers on separate lines -----
getNumbers=function(start, parFile, numSamples) {  
  for (i in 1:numSamples) {
    curnum=as.numeric(unlist(strsplit(parFile[start+i], split=separator))[1])
    if (i==1) {
      num=curnum 
    } else {
      num=c(num, curnum)
    }
  }
  num
}

#--- Get population sizes -----------------
start=3
popSizes=getNumbers(start, parFile, numSamples)
#Rescaling pop sizes
popSizes=round(popSizes*rescalingFactor, digits=0)

iniPopSizes=popSizes

#--- Get sample sizes -----------------

readSampleSizesTimesAndInbreedingLevel=function(start, parFile, numSamples) {
  for (i in 1:numSamples) {
    curLine=unlist(strsplit(parFile[start+i], split=separator))
    curSampSize=as.numeric(curLine[1])
    curSampTime=0
    curInbreeding=0
    if (length(curLine)>1) curSampTime=as.numeric(curLine[2])
    if (length(curLine)>2) curInbreeding=as.numeric(curLine[3])
    if (i==1) {
      sampSize=curSampSize
      sampTime=curSampTime
      inbreeding=curInbreeding
    } else {
      sampSize=c(sampSize,curSampSize)
      sampTime=c(sampTime,curSampTime)
      inbreeding=c(inbreeding,curInbreeding)
    }
  }
  list(ss=sampSize, st=sampTime, inb=inbreeding)
}

start=start+numSamples+1
# sampSizes=getNumbers(start, parFile, numSamples)
sampSizesStats=readSampleSizesTimesAndInbreedingLevel(start, parFile, numSamples)
#Rescaling sample times
sampSizesStats$st=round(sampSizesStats$st*rescalingFactor, digits=0)

sampSizes=sampSizesStats$ss
sampTimes=sampSizesStats$st
inbrCoeff=sampSizesStats$inb

#--- Get growth rates -----------------
start=start+numSamples+1
growthRatesInitial=getNumbers(start, parFile, numSamples)
# save this into growthRates which will be used and updated when printing historical events
growthRates=growthRatesInitial

#--- Get number of migration matrices -----------------
start=start+numSamples+1
numMigMat=as.numeric(unlist(strsplit(parFile[start+1], split=separator))[1])

#--- Read migration matrix
readMigMat=function(start, parFile, numSamples) {
  for (i in 1:numSamples) {
    curmigs=as.numeric(unlist(strsplit(parFile[start+i], split=separator)))
    if (i==1) {
      migs=curmigs 
    } else {
      migs=rbind(migs, curmigs)
    }
  }
  rownames(migs)=1:numSamples
  migs 
}

#--- Get migration matrices as a list --------------
start=start+2
migMats=list()
if (numMigMat>0) {
  for (i in 1:numMigMat) {  
    curMigMat=readMigMat(start, parFile, numSamples) 
    migMats[[i]]=curMigMat
    start=start+numSamples+1
  }
}

#Rescaling migration rates
if (numMigMat>0) {
  for (i in 1:numMigMat) {
    migMats[[i]]=migMats[[i]]/rescalingFactor
  }
}

#--- Get number of historical events
start=start+1
numHistEvents=as.numeric(unlist(strsplit(parFile[start], split=separator))[1])

###################### HISTORICAL EVENTS HANDLING ##############################

#..... Read Historical Event .......
last.he.time=0
if (numHistEvents>0) {
  for (i in 1:numHistEvents) {
    start=start+1
    #Take care of nomig keyword
    nomig=F
    if (grepl("nomig", parFile[start])) {
      nomig=T
      gsub("nomig", "", parFile[start])
    }
    curHE=as.numeric(unlist(strsplit(parFile[start], split=separator)))
    if (nomig) curHE[7]=-1 
    #Rescaling time of event
    curHE[1]=round(curHE[1]*rescalingFactor, digits=0)
    if (i==1) {
      histEvents=curHE
      last.he.time=curHE[1]
    } else {
      histEvents=rbind(histEvents, curHE)
      if (histEvents[i,1] > last.he.time) last.he.time=histEvents[i,1] 
    }
  }
  if (numHistEvents>1) rownames(histEvents)=1:numHistEvents
}
names(last.he.time)=""
last.he.time=as.numeric(last.he.time)

yTimeLimit=0
if (last.he.time!=0) yTimeLimit=last.he.time*(1+propLastsegment) #Add propLastsegment to y axis after last event (to draw stuff)

#Reorder events by their times
if (numHistEvents>1) histEvents=histEvents[order(histEvents[,1],decreasing=FALSE),] else {
  if (numHistEvents==1) histEvents=matrix(histEvents, nrow=1,  byrow = T)
}

endReadParFile=start



##############################  PLOTTING THE MODELED SCENARIO ######################################


#--- Graphical functions ...........................................................................

fullHeadArrow=function(x0, y0, x1, y1, length, angle, color="black",  weight=1) {
  arrows(x0, y0, x1, y1, length, angle, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.80, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.60, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.40, code=2, lty=1, col=color, lwd=weight) 
  arrows(x0, y0, x1, y1, length, angle*0.20, code=2, lty=1, col=color, lwd=weight)
  arrows(x0, y0, x1, y1, length, angle*0.10, code=2, lty=1, col=color, lwd=weight)
}
drawTriangle=function(growth, x, y, size, aspRatio, color) {
  if (growth>0) {
    x0=x; y0=y
    x1=x-size/2; y1=y+size/2*aspRatio
    x2=x+size/2; y2=y+size/2*aspRatio
  } else {
    x0=x-size/2; y0=y
    x1=x; y1=y+size/2*aspRatio
    x2=x+size/2; y2=y
  }
  polygon(c(x0, x1, x2, x0), c(y0, y1, y2, y0), col=color)
  return(y+size/2*aspRatio)
}


#--- Computing maximum current pop size ............................................................
maxPopSize=popSizes[1]
minPopSize=popSizes[1]
if (length(popSizes)>1){
  for (i in 2:length(popSizes)) {
    if (popSizes[i]>maxPopSize) maxPopSize=popSizes[i];
    if (popSizes[i]<minPopSize) minPopSize=popSizes[i];
  }
}

#--- Find min and max pop sizes over the whole population history ..................................
ps=popSizes
isGrowth=FALSE
if (numHistEvents) {
  #Need to keep track ofgrowth rates over time
  gRates=growthRates
  prevTime=0
  for (i in 1:numHistEvents) {
    he=histEvents[i,]
    curTime=he[1]
    sink=he[3]+1 
    resize=he[5]
    growth=he[6]
    
    #Begin by resizing pop sizes due to growth since last event
    for (j in 1:length(ps)) {
      ps[j]=ps[j]*exp(gRates[j]*(curTime-prevTime))
      if (ps[j]>maxPopSize) maxPopSize=ps[j]
      if (ps[j]<minPopSize) minPopSize=ps[j]
      if (ps[j]==0) print(paste("Deme ", j-1, " reaches size zero at time ",  curTime, " due to large negative growth (", gRates[j], ")", sep=""))
      if (is.infinite(ps[j])) print(paste("Deme ", j-1, " reaches infinite size at time ",  curTime, " due to large positive growth (", gRates[j], ")", sep=""))
    }
    prevTime=curTime
    
    #Handle transformed keep keyword
    if(growth!=-9999) {
      gRates[sink]=growth 
    }
    
    #Implement resize
    ps[sink]=ps[sink]*resize
    if (ps[sink]>maxPopSize) maxPopSize=ps[sink]
    if (ps[sink]<minPopSize) minPopSize=ps[sink]
  }
}

#-- Function to compute the circle radius for a given pop size ....................................
interpolRadius=function(curSize, minSize, maxSize, minRadius, maxRadius, logScale) {
  if(logScale) {
    minSize=log10(minSize)
    maxSize=log10(maxSize)
    curSize=log10(curSize)
  }
  curRadius=minRadius+(curSize-minSize)*(maxRadius-minRadius)/(maxSize-minSize)
  curRadius
}

#==============================                     ========================
#==============================   BEGIN MAIN PLOT   ========================
#==============================                     ========================

library("plotrix")
library("diagram")

if (printPDF) {
  pdf(pdfFileName, width=pdf.x.size, height=pdf.y.size)
} 

par(xpd=F, mar=c(6,4,3,0.5))

maxRadius=maxRadius*(numSamples+2)
minRadius=minRadius*(numSamples+2)

title=parFileName

plot(x=1:numSamples, type="n", xlab="", ylab="time (gen)", xlim=c(-0.5, numSamples+1.5), cex.main=0.8,
     ylim=c(0, yTimeLimit), main=title, xaxt = 'n', log=logScaleAxis, cex.axis=0.9, cex.lab=0.9)
axis(side=1, labels=c("Mig Mat",0:(numSamples-1), " \nTimes"), at=0:(numSamples+1), cex.axis=0.8)
mtext("demes", side=1,line=2, cex=0.8)

w <- par("pin")[1]/diff(par("usr")[1:2])
h <- par("pin")[2]/diff(par("usr")[3:4])
aspRatio <- w/h

# aspRatio=yTimeLimit/(numSamples+2)

#--- Plotting population circles according to their current size and sampling time
for (i in 1:numSamples) {
  curRadius=interpolRadius(popSizes[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
  #   print(curRadius)
  if (i==1) topOfCircle=sampTimes[i]+curRadius*aspRatio else topOfCircle=c(topOfCircle, sampTimes[i]+curRadius*aspRatio)
  curColor=popCol
  if (sampSizes[i]==0) curColor="white"
  if (inbrCoeff[i]==0) {
    draw.circle(i, sampTimes[i], radius = curRadius, col=curColor, border=popBorderCol) 
  } else {
    floating.pie(i, sampTimes[i], c(inbrCoeff[i], 1-inbrCoeff[i]), radius = curRadius, 
                 col=c(inbreedColor,curColor), startpos=pi, border=popBorderCol)    
  }
  if (sampTimes[i]>0) {
    text(i, sampTimes[i]-curRadius*aspRatio, labels=sampTimes[i], cex=timeProp, col=oldSampColor, pos=1)
  }
  #Draw a vertical arrow in case of pop growth
  curGrowthRate=growthRates[i]
  
  #--- Draw gtrowing or shrinking triangle on top of pop circle to reflect growth type
  if (curGrowthRate!=0) {
    arLength=0.15
    topOfCircle[i]=drawTriangle(curGrowthRate, x=i, y=topOfCircle[i], arLength, aspRatio, color=growthCol)
  }
}

#--- Handle first migration matrix .......................
curMigMatNum=0
curvature=0.0075*last.he.time
text(0-migrOffset, 0, labels=0, cex=migMatNameProp, col=migrMatCol)
if (numMigMat) {
  curMigMat=migMats[1][[1]]
  for (sink in 1:numSamples) {
    for (sourc in 1:numSamples) {
      if (sink!=sourc & curMigMat[sourc, sink]>0)  {
        differ=sourc-sink
        curvedarrow(from=c(sourc, 0), to=c(sink,0), curve=curvature*(abs(differ)*0.55^abs(differ)), arr.adj=1, 
                    arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY, 
                    lcol=migrMatCol, arr.length=arrowLength)   
        if (plotMigrRates) {
          curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)        
          if (differ>0) {
            text(sink+abs(differ)/2, aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
          } else {
            text(sourc+abs(differ)/2, -aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
          } 
        }
      }
    }
  }
}

#---Draw all events on the population tree
lastTime=0
activePops=1:numSamples
numActivePops=numSamples
lastSink=-1
if (numHistEvents) {
  for (i in 1:numHistEvents) {
    #Extract historical event
    he=histEvents[i,]
    he.time=he[1]
    he.source=he[2]+1 #+1 is due to the use of base 0 for deme number in fsc
    he.sink=he[3]+1
    he.migr=he[4]
    he.resize=he[5]
    he.growth=he[6]
    if(he.growth==-9999) { #handle transformed keep keyword
      he.growth=growthRates[he.sink] #Keep current growth rate
    }
    he.migrMat=he[7]
    
    if(he.migrMat==-9999) { #handle transformed keep keyword
      he.migrMat=curMigMatNum
    }
    
    #-- Draw event time
    if (i%%2) {
      slide=timeOffset
    }
    else {
      slide=-timeOffset
    }
    
    #Draw all vertical segments ....................
    for (p in 1:length(activePops)) {
      #Handle population with older sampling times
      if (sampTimes[activePops[p]]>0) {
        minTime=max(lastTime, topOfCircle[activePops[p]]);
      } else {
        minTime=lastTime;
      }
      if (he.time>sampTimes[activePops[p]] & he.time>topOfCircle[activePops[p]]) {
        if (activePops[p]==lastSink | i==1) {
          segments(activePops[p], topOfCircle[activePops[p]], activePops[p], he.time)    
        } else {
          segments(activePops[p], minTime, activePops[p], he.time)        
        }
      }
    }     
    
    #Handle growth rate changes since last event ......................
    #Update pop sizes according to current growth rates
    if (i>1) {
      prev.he=histEvents[i-1,] 
      branchLength=he.time-lastTime
      for (p in 1:length(activePops)) {
        curPop=activePops[p]
        popSizes[curPop]=popSizes[curPop]*exp(growthRates[curPop]*branchLength);
      }
    }
    #Update growth rate
    growthRates[he.sink]=he.growth
    
    
    lastTime=he.time
    #Handle resize of sink pop ........................
    popSizes[he.sink]=popSizes[he.sink]*he.resize
    
    curRadius=interpolRadius(popSizes[he.sink], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
    #     popRadius[he.sink]=curRadius
    topOfCircle[he.sink]=he.time+curRadius*aspRatio
    draw.circle(he.sink, he.time, radius = curRadius, col=popCol, border=popBorderCol)
    
    
    #--- Draw growing or shrinking triangle on top of pop circle to reflect growth type
    if (he.growth!=0) {
      arLength=0.15
      topOfCircle[he.sink]=drawTriangle(he.growth, x=he.sink, y=topOfCircle[he.sink], arLength, aspRatio, color=growthCol)
    }
    
    #Handle population fusion .........................
    if (he.migr>=1 & he.sink!=he.source) { #This is a population fusion
      if (numActivePops==numSamples) removedPops=he.source else removedPops=c(removedPops, he.source)
      numActivePops=numActivePops-1
      activePops=(1:numSamples)[-removedPops]
      #Draw connecting arrows from source to sink
      fullHeadArrow(he.source, he.time, he.sink, he.time, length=0.15, angle=20)      
      #Redraw time with the right color
      text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=popFusionColor)
    } else {
      #--- Handle admixture event ........................
      if (he.migr>0 & he.migr<1) { 
        #Draw connecting arrows from source to sink      
        segments(he.source, he.time, he.sink, he.time, col=admixCol, lty=2)
        if (he.sink>he.source) {
          fullHeadArrow(he.sink-0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
        } else {
          fullHeadArrow(he.sink+0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
        }     
        #Redraw time with the right color
        text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=admixCol)
      }
      else text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=timeCol)
    }    
    
    
    #--- Handle migmat change  ............................
    if (i!=numHistEvents) nextTime=histEvents[i+1,1] else {
      nextTime=yTimeLimit
    }
    time2DrawArrows=(he.time+nextTime)/2
    
    if (he.migrMat!=curMigMatNum) {
      if (he.migrMat>-1) {
        curMigMat=migMats[he.migrMat+1][[1]]
        for (sink in 1:numSamples) {
          for (sourc in 1:numSamples) {
            if (sink!=sourc & curMigMat[sourc, sink]>0)  {  
              differ=sourc-sink
              curvedarrow(from=c(sourc, time2DrawArrows), to=c(sink, time2DrawArrows), curve=curvature*(abs(differ)*0.55^abs(differ)), 
                          arr.adj=1, arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY, 
                          lcol=migrMatCol, arr.length=arrowLength)  
              if (plotMigrRates) {
                #Write Nm values
                curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)
                if (differ>0) {
                  text(sink+abs(differ)/2, time2DrawArrows+aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
                } else {
                  text(sourc+abs(differ)/2, time2DrawArrows-aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol) 
                } 
              }
            }
          }
          # curMigMatNum=he.migrMat
        }
      } 
      #--- Draw separation between migration matrices numbers on the left
      segments(-migMatLineLength/2, he.time, migMatLineLength/2, he.time, lty=3, col=migrMatCol)
    } 
    curMigMatNum=he.migrMat
    
    #Output current valid migration matrix      
    if (i%%2) {
      slide=migrOffset
    }
    else {
      slide=-migrOffset
    }   
    if (he.migrMat>-1) {
      migText=he.migrMat
      curCex=migMatNameProp
    } else {
      migText="nomig"
      curCex=migMatNameProp/2
      slide=slide*2
    }
    text(0+slide, time2DrawArrows, labels=migText, cex=curCex, col=migrMatCol)    
    
    lastSink=he.sink
  }
}

#-- Draw last branch
segments(activePops[1],topOfCircle[activePops[1]], activePops[1], yTimeLimit)

#==============================   PLOT LEGENDS IN MARGINS   ========================

#Compute space available in margin
minY.coo=grconvertY(0, from="nic", to="user")

par(xpd=NA)
#--- Draw population size scale with circles of different sizes .................
maxOrder=ceiling(log10(maxPopSize))
minOrder=floor(log10(minPopSize))
popSizeRadius=10^(maxOrder:minOrder)
winWidth=numSamples+2
ypos=3/4*minY.coo
text(x=-winWidth/10*1.2, y=ypos, labels="Pop. \nsizes ", cex=.8, pos=2)

for (i in 1:length(popSizeRadius)) {
  curRadius=interpolRadius(popSizeRadius[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
  #   print(curRadius)
  if (curRadius>0) {
    xpos=-winWidth/10+(i-1)*winWidth/10
    draw.circle(xpos, ypos, radius=curRadius, col=popCol, border=popBorderCol)
  }
  text(xpos, ypos-abs(ypos)*0.1, popSizeRadius[i], cex=0.7, pos=1, col="black")
}

#--- Legend for growing or shrinking populations ...............................
if (isGrowth) {
  x=winWidth-1.5*winWidth/10; y=ypos+abs(ypos)*0.2
  text(x, y-abs(ypos)*0.1, labels="Populations", cex=0.8)
  x=winWidth-2*winWidth/10; y=ypos-abs(ypos)*0.1
  drawTriangle(1, x, y, size=0.15, aspRatio, color=growthCol)
  text(x, y-abs(ypos)*0.1, labels="growing", pos=NULL, cex=0.7)
  x=winWidth-winWidth/10
  drawTriangle(-1, x, y, size=0.15, aspRatio, color=growthCol)
  text(x, y-abs(ypos)*0.1, labels="shrinking", pos=NULL, cex=0.7)  
}


#if (printPDF) dev.off()
dev.off()

## null device 
##           1